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Abstract 


We  consider  the  task  of  numerically  approximating  the  solution  of  an  ordinary 
differential  equation  initial-value  problem.  We-art,  interested  in  two  questions: 

(i.)  For  any  given  one-step  method,  what  is  the  complexity  of  finding  an 
approximate  solution  with  error  less  than « ? 

(2.)  Given  an  infinite  sequence  of  one-step  methods  of  increasing  order,  how 
should  the  method  and  the  step-size  be  picked  so  as  to  minimize  the 
complexity  of  finding  such  an  approximation? 

We-  describe  a methodology  that  handles  both  questions.  Furthermore,  we  find  that 

within  such  a sequence  of  methods,  the  following  hold  under  very  general 

circumstance^, 

(1.)  For  any  «,  0 < • < 1,  there  is  a unique  choice  of  order  and  step-size  which 
minimizes  the  complexity. 

(2.)  As  t decreases,  both  the  optimal  order  and  the  complexity  increase 
monotonically,  tending  to  infinity  as  • tends  to  zero. 

These  ntsults  are  applied  to  several  classes  of  one -step  methods.  In  doing  so,  we 

exhibit  some  new  Taylor  series  methods  that  are  asymptotically  bett'tr  than  Runge- 

Kutta  methods  for  problems  of  small  dimension.  Moreover,  we  provu  that  among  all 

classes  of  nonlinear  Runge-Kutta  methods,  those  due  to  Brent  have  the  highest  order 

possible. 
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With  few  exceptions,  past  work  in  analytic  computational  complexity  has  focused 
on  the  problam  of  finding  a zaro  of  a (nonlinear)  transformation  of  Banach  spaces:  in 
most  work,  this  problem  is  specialized  to  that  of  finding  a zero  of  an  operator  on  a 
finite-dimensional  real  or  complex  vector  space  (and  in  much  of  this  work,  the  problem 
is  further  specialized  to  the  one -dimensional  case).  Much  has  been  discovered  about 
the  computational  aspects  of  iterative  schemas  for  the  solution  of  such  problems, 
especially  in  the  areas  of  minimal  complexity  (•.&,  Ktmg  and  Treub  [73i  Traub  and 
Wofniakowski  [76])  and  maximal  order  (e.ju  Kung  and  Traub  [74J  Woiniskowsfci  [75]}, 

In  this  iper,  we  will  consider  another  topic  in  analytic  computational  cowpJexity 
theory,  that  of  finding  complexity  bounds  for  the  numerical  solution  of  ordinary 
differential  equation  initial-value  problems  on  a fixed  Interval.  We  wilt  not  be 
interested  in  questions  of  the  existence  and  the  uniqueness  of  the  solutions  to  such 
problems:  in  fact,  we  will  restrict  our  discussion  of  the  application  of  ganerel  results  to 
the  case  where  the  unique  solutions  to  these  problems  are  analytic  functions. 

We  will  limit  oursalves  here  to  classes  of  one-step  methods  tor  the  numeric  &1 
solution  of  these  problems;  in  terms  of  informational  usage,  the&a  methods  are 
analogous  to  iterative  zero-finding  methods  without  memory  (Traub  [64J  [72]K 
Analogous  to  the  one-point  itarative  methods  with  memory  are  the  wuiHsteo  methods 
for  initial-value  problems;  these  methods  will  be  dealt  with  in  $ future  paper. 


Our  approach  will  be  'o  assume  that  an  initial-value  problem  is  given,  along  with 
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soma  error  criterion  i,  where  0 < c < lj  we  then  wish  to  compute  an  approximate 
solution  with  error  no  greater  than  «.  Two  basic  questions  concern  us: 

(1.)  For  any  given  method,  what  is  the  complexity  of  solving  this  problem? 

(2.)  Given  any  "basic*  sequence  of  methods  with  increasing  order,  which  method 
has  minims!  con  plexity? 

In  Section  2,  we  describe  a methodology  that  handies  both  questions  for  classes 
consisting  of  methods  whose  error  functions  have  a special  form.  Furthermore,  we 
find  that  within  such  a basic  sequence  of  methods,  the  following  hoSd  under  very 
general  conditions: 

(1.)  For  any  i,  there  is  s unique  choice  of  order  and  step  size  minimizing  the 
complexity. 

(2.)  As  « decreases,  both  the  optima]  order  and  the  complexity  increase 
monotorocaiiy,  tending  to  infinity  as  a te. *d$  to  zaro. 

Furthermore,  within  r«*ny  classes  of  problems  and  methods,  the  "penalty"  (o4F»  the 

amount  th®  cost  curve  turns  r»ar  the  optimum)  associated  with  using  non-optima! 

order  tends  to  infinity  as  s tends  to  zero. 

These  conclusions  are  an  interesting  contract  to  known  resuits  on  zero-finding 
via  iterations  without  memory.  The  latter  results  tend  to  support  the  "folklore"  idea 
that  it  is  "bettor"  to  use  a low-order  method  msny  times,  than  to  use  a high-order 
method  a few  times.  In  the  -paint  case,  ootimai  order  is  low,  while  in  the  multipoint 
case,  optimal  &rdar  increases  with  the  problem  complexity  (but  with  little  penalty  for 
using  a method  of  non-optimal  ortfer)  (Kung  and  Trsub  [73]).  In  addition,  optimal  order 
for  these  pcofeism  does  not  depend  on  the  error  criterion!  it  is  computed  tar  the 
limiting  case  as  « approaches  z&ro. 

One  may  wonder  why  there  is  this  discrepancy  between  the  resultc  for  the 
initial-value  problem  and  thosa  lor  tbs  zoro-finding  problem,  since  any  initial-value 
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problem  may  bo  written  si  an  operator  equation,  as  in  Stetter  [73}  The  reason  for 
this  is  that  the  methods  usad  for  the  two  problems  differ  greatly— those  for  the  iratial- 
vaiue  problem  compute  estimates  for  solution  values  at  new  points  by  discretization. 
while  those  for  zero-finding  compute  improved  estimates  for  the  zero  of  a function  by 
iteration- 

In  Section  3,  we  discuss  the  extension  of  these  results  to  daises  consisting  of 
methods  whose  error  functions  are  somewhat  more  complicated  than  those  considered 
in  Section  2. 

In  Section  4,  we  introduce  the  nations  of  normality  and  order -convergence  for  a 
basic  sequence  of  one-step  methods.  We  prove  that  they  are  equivalent  under  certain 
circumstances.  A basic  sequence  of  methods  enjoying  these  properties  is  very  easy  to 
deal  with  in  many  respects,  especially  whan  one  is  interested  in  comparing  upper  and 
lower  complexity  bounds  for  such  a class. 

In  Sc  ‘ion  5,  we  apply  the  theory  developed  in  th«  proceeding  sections  e the 
general  problem  of  an  autonomous  system  of  equations.  We  show  that  the  optimal 
order  and  complexity  behave  as  described  by  (I.)  and  (2.)  above  for  the  class  of 
Taylor  series  methods  and  for  various  classes  of  Runge-Kutta  methods.  In  addition,  we 
construct  new  Taylor  series  methods  that  are  asymptotically  bet  tar  (as  « tends  to 
zero)  than  Runge-Kutta  methods  for  problems  of  small  (^mansion. 

In  Section  6,  we  look  at  the  problem  of  a single  scalar  autonomous  equation.  In 
this  case,  we  may  use  the  classes  of  "nonlinear  Runge-Kutta  methods"  developed  by 
Brent  [74}  [76}  We  show  that  the  behavior  described  by  (1.)  and  (2.)  above  holds  for 
these  methods.  In  addition,  we  prove  that  among  all  classes  of  nonlinear  Runge-Kutta 
methods,  those  due  to  Brent  have  the  highest  order  possible. 


Section  7 describes  some  numerical  data  that  support  the  above  theoretical 
results.  In  particular,  these  data  seem  to  indicate  that  even  for  modest  values  of  t, 
there  are  considerable  savings  in  using  methods  of  optima!.,  rather  than  fixed,  order. 

Finally,  in  Section  8,  we  draw  some  conclusions,  make  some  comparisons,  point 
out  some  unanswered  questions,  and  define  new  areas  to  which  this  theory  should  be 
extended. 


Section  2 


Optimality  Within  a Strong  Basic  Sequence 

We  are  interested  in  the  numerical  solution  of  a class  of  ordinary  differential 
equation  initial-value  problems  on  a fixed  interval  I of  finite  length}  we  tako  I - [0,  1] 
without  loss  of  generality.  More  precisely,  let  5)  be  a set  of  initial-value  points  in  the 
real  N-dimensional  linear  space  R^ , and  let  be  a set  of  operators  on  R^ , such  that 
the  Initial-value  problem  of  finding  a function  x:  I •*  R^  satisfying 

x(t)  - v(x(t))  if  t € int  I , 

(2.1) 

x(0)  - x0 

has  a unique  solution  for  every  (xa  , v)  « The  autonomous  form  of  this  system 

is  no  restriction,  since  any  ncn-autonomous  system  may  be  made  autonomous  by 
increasing  the  dimension  of  the  system  by  one.  • 

The  moaei  of  computation  to  be  used  is  fairly  general.  We  assume  only  that  all 
arithmetic  operations  are  performed  exactly  in  R (f.e.,  infinite-precision  arithmetic), 
and  that  for  any  algorithm  to  be  considered  for  the  solution  of  (2.1),  a set  of 
procedures  is  given  for  the  computation  of  any  information  about  v squired  by  that 
algorithm.  (For  instance,  with  Runge-Kutta  methods,  we  must  be  able  to  compute  v at 
any  point  in  its  domain.) 

In  this  paper,  we  are  interested  in  the  numerical  solution  of  (2.1)  via  one-step 
methods,  using  an  equidistant  grid  as  defined  in  Stetter  [73].  (We  limit  ourselves  to 
equidistant  grids  in  order  to  facilitate  the  comparison  of  methods  of  different  orders; 
the  other  extreme  is  taken  by  Lindberg  [74],  who  considers  the  problem  of  picking  an 
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optimal  grid  for  a given  method  of  fixed  order.)  Thus  the  methods  considered  will 
generate  approximations  Xj  to  x(t|)  by  the  recursion 

(2.2)  X|+i  - xj  + h ^Xj  ,h)  (Osi£n-l,n-  h“*) 

where  h is  the  step-size  and  f>  is  the  increment  function  for  the  method  (Hanrici  [62]h 
for  briefness,  we  will  refer  to  "the  method  f."  Despite  the  fact  that  ?(xj , h)  will 
depend  on  some  information  about  v,  we  will  not  explicitly  indicate  this  dependence. 

Thus,  the  method  ? produces  an  approximation  to  the  true  solution  of  (2.1).  We 
want  to  measure  the  discrepancy  between  the  approximate  and  true  solutions.  Various 
error  measures  have  been  introduced  in  the  literature.  These  include  the  local 
truncation  error  per  step,  the  local  truncation  error  per  unit  step,  and  the  global  errors 
see  Henrici  [62]  or  Stetter  [73]  for  definitions.  These  error  measures  may  be  either 
absolute  or  relative  (in  the  usual  sense);  they  may  be  measured  either  at  the  endpoint 
of  the  interval  (as  in  Henrici  [621  Hindmarsh  [74])  or  over  the  entire  grid  (a3  in 
Sandberg  [671  lindberg  [74]).  There  has  been  a great  deal  of  discussion  of  which 
error  criterion  is  the  best  one  to  use;  for  instance,  Gear  [71]  (Section  9.3)  uses  local 
error  per  step,  white  Hull  et  al.  [72]  use  local  error  per  unit  step.  We  take  no  sides  in 
this  discussion,  since  any  of  these  error  measures  may  be  used  in  the  analysis  to 
follow. 


Before  proceeding  any  further,  we  will  establish  some  notational  conventions. 
Let  DC  be  an  ordered  ring;  then  X * and  DC++  will  respectively  denote  the  nonnegative 
and  positive  elements  of  DEL  (This  will  be  used  in  the  cases  DC  ••  R,  the  real  numbers, 
and  DC  - Z,  the  integers.)  The  symbol  means  "is  defined  to  be,"  while  V means 
"is  identically  equal  to."  The  symbol  "V"  will  be  used  to  _<snote  the  gradient  of  a 
.mapping.  If  xi»  X2:  **  **  R and  R^  -»  R are  differentiable,  then  for  i ■ 1,  'i,  we  will 
write 
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dj  «^XlO),X2^)) 

for  the  result  of  differentiating  w'^  respect  to  Xj»  and  than  substituting 

X^  - Xi<t),  X2  “ 'X2^  us®  notations  "x  i a"  and  "x  t a"  to  indicate  one-sided 
limits  as  in  Bu' « [65].  Finely,  we  shall  write  "(a.b)c"  to  inidicate  the  c*h  part  of 
equation  (a.b),  as  in  Gurtin  f/5]. 

Now  we  are  prepared  to  define  our  problem.  Let  2)  and  7?  bs  as  above;  ^ 
consider  a problem  (xq,v)  in  JDx^.  Let  * be  a class  of  one-step  methods,  and  let 
9:  4x1  -*  R+  satisfying  lim  e{?,h)  - 0 be  a given  function  that  will  serve  as  an 
error  measure.  Choose  an  error  criterion  t satisfying  the  technical  restriction  0 < « < 

1.  We  then  wish  to  answar  two  questions: 

(1.)  Given  * < 4,  how  may  we  pick  h < I such  that 
(2,3)  #<*h)  5 ., 

and  what  is  the  complexity  of  the  process  defined  by  p and  h? 

(2.)  How  may  one  choose  among  all  (*,h)  « 4XI  such  that  (2J3)  holds,  that  pair 
(**,h*)  giving  minimal  complexity? 

In  order  to  get  useful  bounds  on  *(*,b),  it  is  necessary  to  introduce  the  concept 
of  order.  In  this  section,  we  will  use  a highly  restricted  definition,  which  we  wili  relax 
in  Section  3.  Let  4 - (yp:  p € Z++},  and  suppose  that  there  is  an  analytic  function 
R+  -*  R+  such  that  lim  p_^j  «(p)^P  exists  and  is  nonzero  and 

(2.4)  e(*p,h)  - (rip)  hP  for  h < I and  p < Z++  . 

Then  ^p  is  said  to  have  strong  order  p with  respect  to  *,  and  4 is  said  to  be  a strong 


basic 


j.  (Although  the  error  coefficient  * will  generally  depend  on  the  solution 


x of  (2.1),  we  do  not  explicitly  indicate  this  dspendence.)  Note  that  the  order  of  a 
method  depends  on  the  error  measure;  for  example,  the  order  with  respect  to  the  local 
error  per  step  is  one  greater  than  that  with  respect  to  the  local  error  per  unit  step  or 
the  global  error. 
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Equation  (2.4)  is  somewhat  more  restrictive*  than  that  which  is  usually 
encountered  in  practice;  more  often,  we  expect  « to  depend  on  h.  We  consider  the 
extension  of  our  results  to  this  case  in  the  next  section. 

We  now  are  able  to  measure  the  complexity  of  computing  an  approximate 
solution  to  (2.1),  with  error  not  exceeding  «,  using  a strong  basic  sequence  ♦.  Indeed, 
(2.4)  implies  that  a necessary  and  sufficient  condition  for  «p{f»p,h)  - s is  that 


where 


h - h(p^r)  «(p)"*/P  e“*/p  , 


« In  (i-^) 


(Note  that  since  0 < « < 1,  we  have  « € R**.)  Thus,  the  number  of  steps  needed  is 


given  by 


n if*  - «Kp)^p  e*^p  . 


(Note  that  n (as  given  by  (2.7))  need  not  be  an  integer.  But  this  poses  no  essential 
difficulty;  see  (e.g.)  Traub  and  Woiniakowski  [76})  Next,  suppose  that  there  exists  an 
analytic  function  c;  R + -♦  fl  + such  that  c(p)  is  the  cost  per  step  associated  with  the 
method  *p.  Finally,  we  assume  that  the  cost  per  step  does  not  vary  from  step  to  step; 
for  the  classes  of  methods  we  consider,  this  means  only  that  we  assume  that  the  cost 
of  evaluating  v (or  its  derivatives)  does  not  depend  on  the  point  of  evaluation.  Thus 


given  by 


<ity  C(p/t)  of  solving  (2.1)  to  within  an  error  criterion  « - e"*  is  simply 


C(p^t)  - n c(p)  - f(p)  e*/P  , 


where  we  define  f;  R + -*  R+  by 

(2.9)  f(p)  dp)1^  c(p)  . 

We  now  turn  to  the  question  of  picking  for  each  a € R**  that  order  p giving 
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minimal  complexity.  In  the  analysis  to  follow,  we  will  drop  the  restriction  that  p must 
be  an  integer.  However,  we  will  recover  optimality  over  the  integers  from  optimality 
over  the  real  numbers  in  Corollary  2.1  . Without  loss  of  generality,  we  assume  that 


(2.10) 


p > 0 implies  f(p)  > 0 . 


(If  there  were  a p > 0 with  f(p)  « 0,  use  of  the  method  ?p  would  yield  a solution  with 
zero  complexity,  i.e.,  "with  no  effort.")  In  addition,  wu  assume  that 


(2.11) 


liWpro,  >'(p)  - 


By  (2.9),  this  assumption  maybe  viewed  as  a simple  consequence  of  two  conditions, 
both  of  which  are  quite  natural.  The  first  is  that  iimp^g,  c(p)  - +o;  the  "better"  a 
method  is  (i.e.,  the  higher  its  order  is),  the  more  we  should  expect  to  pay  for  its  use. 
The  second  condition  is  that  if  iim  p^  «(p)  - 0 , then  there  must  exist  a ft  € 1 surh  that 
<t(p)  i for  p sufficiently  large.  (For  example,  in  the  class  of  Taylor  series  methods, 
using  the  worst-case  local  error  per  unit  step  as  the  error  measure,  this  second 
condition  would  follow  from  the  assumption  that  any  problem  (xq,v)  € must  have 
an  analytic  solution.) 

Thus  in  order  to  find  a minimum  for  C(  * ,«),  we  merely  differentiate  (2.8)  with 


respect  to  p,  finding 
(2.12) 


C(p,«)  - p-2  f(p)  e*/P  [G(p>  - «]  , 


where  G:  R ++  -»  R is  given  by 

(2.13)  G(p)  p2  f'(p)/f(P>  . 

Thus  a necessary  condition  that  p be  a minimum  for  C(  • ,*)  is  that  dj  C(p^r)  - 0,  i.e., 


(2.14) 


G(p)  ■ a . 


Sufficient  conditions  for  the  existence  and  uniqueness  of  a p satisfying  (2.14)  and 
minimizing  C(  • ,a ) are  given  in 
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Theorem  2.1;  Let  f satisfy  (2.10)  and  (2.11).  Suppose  also  that 


(2.15) 


G'(p)  > 0 whenever  G(p)  > 0. 


Then  there  is  a function  p*:  I?**  -*  F44,  such  that  (2.14)  holds  if  and  only  if  p - p*(«). 


Moreover,  for  all  p € R -M‘, 
(2.16) 


C*(«)  C(p*(a),or)  S C(p,«), 


with  equality  holding  if  and  only  if  p - p*(«). 

(Since  p*(«)  satisfies  (2.15),  we  call  p*(«)  the  optimal  order.  C*(«)  the  optimal 


(2.17) 

the  optimal  step-size.) 


h*(«)  h(p*(a),ot) 


Proof  of  Theorem  2.1:  If  we  write  the  Maclaurin  series  of  f and  substitute  it 


into  (2.13),  it  is  easy  to  see  that 
(2.18) 


limp+0G(p)  - 0. 


We  now  claim  that 


(2.19)  limpToo  G.(p>  - -ko. 

Indeed,  since  (2.11)  holds  there  is  a pq  > 0 such  that  f'(pQ)  > 0,  i.e.,  G(pq)  > 0.  Thus 
by  (2.15),  G is  monotone  increasing  on  [pq,  +co),  and  hence  either  (2.19)  holds  or  there 
exists  a y > 0 such  that  limpid  G(p)  « 7.  If  the  latter  holds,  then  G is  bounded,  and 


we  nave 


f'(t)/f(t)  S It-2  (1  S t<+«o) 


for  some  4 > 0;  integrating  the  above  inequality  over  1 £ t S p yields 

f(p)  S f(l)a6(1  ‘ l/P>, 

so  that  limp^  f(p)  S f(i)  e®,  contradicting  (2.11).  Thus  (2.18)  and  (2.19)  hold; 
together,  they  imply  that  for  any  a > 0,  there  is  a choice  of  p such  that  (2.14)  holds. 
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Suppose  that  for  some  a > 0,  there  were  two  numbors  pq  < pj  with  (Kpg)  ■ 
G(pj)  - a.  Then  by  Rolle’s  Theorem,  there  is  a P2  between  pq  and  pj  with  G'(p2>  • 0, 
contradicting  (2.15).  Thus  for  each  a > 0,  there  is  a unique  choice  of  p such  that  (2.14) 
holds)  we  denote  this  choice  by  p*(«). 

To  prove  (2.16),  differentiate  (2.12)  with  respect  to  p to  find 
(2.20)  dx2  C(p^r)  - p"’  f(p)  e«/P  G'(p)  + [G(p)  - «]  (d/dp)  [p“2  f(p)  e“/P]  . 
But  upon  substituting  p - p*(a),  the  second  term  in  (2.20)  vanishes  and  the  first  term 
is  positive;  so  we  have 

dj2  C(p*(a),a)  > 0 . 

Thus  p*(«)  gives  a local  minimum  for  C(  • ,a),  which  has  only  one  critical  point  (since 
(2.14)  has  a unique  solution)  and  (2.16)  follows.  | 

Note  that  we  have  not  said  that  p*(«)  is  an  integer;  in  fact,  this  need  not  be  true 
in  general.  Since  the  basic  sequence  $ is  indexed  by  Z**,  we  have  not  yet  solved  the 
problem  of  choosing  from  among  all  (*plh)  such  that  (2^)  holds,  that  pair  yielding 
minimal  complexity.  This  problem  is  solved  by 

Corollary  2.1:  For  any  a > 0,  define  p**(«)  € Z+*  be  t^at  element  of  the  set 
{[p*(tf)J  , fp*(«)"|}  which  gives  the  smaller  value  of  C(  • ,«).  The.'* 

C(p**(c),a)  £ C(p^t)  for  p « Z^. 
with  equality  if  and  only  if  p - p**(«). 

Proof:  Clearly  we  need  only  consider  the  case  where  p*(«)  is  not  an  integer. 
Suppose  there  exists  pq  € Z 4"\  not  equal  to  p**(«),  with  C(pg,e)  S C(p**(«),a).  Without 
loss  of  generality,  assume  pg  < jp*(c}J.  Then  C(pgyr)  i C([p*(e)i«)  i.  C(p*(<r),e),  which 
implies  that  there  is  a pj  € (pg  , p*(a))  such  that  dj  J(pjys)  - 0.  Hence,  G(pj)  - a,  but 
pj  i4  p*(«).  This  contradicts  Theorem  2.1.  | 


1 


It  may  be  readily  verified  that  the  hypotheses  2.1  are  satisfied  for 

many  classes  of  functions  f.  Some  of  these  are 
t rythmic:  f(p)  - In  (p  + a) , 
monomial:  f(p)  - pm  (m  i R **) , 
exponential:  f{p)  - fiP  (0  > 1) , 
super-exponential:  f<p)  » pP  , and 
hyper-exponential:  f(p)  - . 

(Wtj  write  "In  (p  + e)",  where  e is  the  base  of  the  natural  logarithms,  rather  than  "In  p" 
as  a technical  convenience.  However,  an  expression  of  the  form  "In  (p  + f)H  with  > 0 
is  necessary  to  guarantee  that  f(l)  > 0.)  Furthermore,  we  find  that  if  f to  the 
monomial-logarithmic  form 

f(p)  - p*  On  -*))b  (a,  b ( R ^) , 

then  the  hypotheses  of  Theorem  2.1  hold.  This  may  be  verified  either  directly,  or  by 
using  the  following  Lemma,  along  with  the  tact  that  the  hypotheses  hold  for  f(p)  « p 
and  f(p)  - In  (p  + e). 

Lemma  2.1:  Let  f hive  the  form 

f(p)  - a n£j  (fjfp)/1  , 

where  a € R**,  and  for  each  i (1  S 5 £ m),  fj  satisfies  the  hypotheses  of  Theorem  2.1 
and  rj  f R++.  Then  f satisfies  the  hypotheses  of  Theorem  2.1. 

Proof:  !t  is  clear  that  if  each  fj  satisfies  (210)  and  (211),  then  so  does  f.  If 
each  fj  yields  (via  (213))  a Gj  satisfying  (215),  the:,  f yields  a G in  the  form 

G(p)  - I?!  rjGj(p) , 

and  so  Vi  is  clear  that  G sattsfies  (215).  2 

Note  that  for  our  purposes,  we  will  only  be  interested  in  monomial  and 
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mwvsomial-logarithmic  growth.  We  include  the  other  examples  of  functions  that  satisfy 
the  hypotheses  of  Theorem  2.1  to  illustrate  the  wide  variety  of  functions  that  qualify. 

So  we  havr>  seen  that  under  the  hypotheses  of  Theorem  2.1,  there  is  a unique 
choice  of  order  ind  step  size  minimizing  ? e total  complexity  i^r  any  error  cr'lerion. 
What  happens  to  these  choices  as  a changes? 

Theorem  2.2:  Let  f satisfy  the  hypotheses  of  Theorem  2.1.  Then 

(1.)  p*(«)  and  C*(«)  increase  monotonically  with  e. 

(2.)  Iima|03  p*(«)  - linigfto  C*(«)  - +oo. 

(3.)  If  there  exists  M > 0 such  that  «(p)^p  £ to  for  all  p,  then 
lim  inf^Qj  h*(«)  > 0 if  e/p *{«)  is  bounded  as  atco. 

Proof:  To  prove  (1.),  note  that  p*  is  the  functional  inverse  of  G.  Thus  p*'(«)  « 
G'{p*(«))_1  > 0,  so  that  p*(«)  increases  with  ...  Now  use  the  chain  rule: 

C*'(«)  - 6lC(p*im)^p*'{m)  <•  ^ C(p*(a),a). 

But  the  first  term  on  the  right-hand  side  vanishes  by  the  definition  nf  p*{«).  So 
C*'<«)  - a2  C(p*(«U)  - (p^e))'1  f(p*(«))  »«/>*<«>  > 0 
ar*J  C*(a)  incre*\i«ss  with  a. 

Suppose  that  lim  p*(a)  / *co  . Sinca  p*(a)  increases  monotonically  with  a, 
there  is  an  L > 0 such  that  lim^gj  p*(a)  - L So  (2.14)  implies  thet 
G(L)  ^ limafTO  G<p*(a))  - lim^too  « “ +*° » 
contradicting  the  continuity  cf  G.  This  proves  the  first  part  of  (2.)  . Now  for  any 
« > 0,  we  have 

C*(a)  - {(p*{«))  e*/P^*>  > f(p*(«)) . 

Let  «T co;  then  (2.1 1)  and  the  first  part  ct  (2.)  imply  that  the  second  part  of  (2.)  holds. 

To  prove  (3.),  let  such  an  to  > 0 exist,  so  that  [h*(«)]“*  £ M e*/p  ^a\  Then  we 
see  that  lim  inf^Q,  h*(a)  > 0 if  [h*(a)]"*  r bounded  as  aTco,  which  itself  is  true  if 
a/p*(a)  is  bounded  as  atco,  | 
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Therefore  under  a vary  general  sat  of  conditions,  ws  sea  that  the  more 
accuracy  we  want  in  our  computed  solution,  J.he  greater  its  complexity  becomes.  Of 
cour53,  this  is  just  what  we  would  expect.  What  is  Kurnev  hat  surprising  is  that  the 
minmal  complexity  is  obtained  by  letting  the  order  ,>  increase  as  the  error  i 
decreases,  with  p increasing  without  bound  as  « tends  to  zero.  Moreover,  the  last  part 
of  the  theorem  says  that  not  only  should  the  ordor  be  increased  when  trying  to  obtain 
a more  accurate  solution,  but  that  it  may  dually  turn  out  that  the  step-size  should 
not  be  allowed  to  tsnd  to  zero. 

We  now  determine  whether  we  are  saving  a great  deal  by  using  the  optimal- 
order  method.  This  may  be  thought  of  in  several  ways;  we  will  consider  how  sharply 
the  cost  curve  turns  at  the  optimum,  the  cost-difference  between  is,.ng  a method  of 
fixed  order  and  a method  of  optimal  order,  and  the  cost-ratio  of  a fixed-order  method 
to  an  optimal-order  method.  We  will  show  that  under  certain  reasonable  conditions,  all 
cf  these  measures  tend  to  infinity  with  m. 

How  sharply  the  cost  curve  turns  at  the  maximum  is  measured  by  C(p*(«),«). 
If  W'o  consider  five  of  the  growth  models  mentioned  above  (e^,  monomial,  monomial- 
logarithmic,  exponential,  hyper-exponential,  and  super-exponential),  we  find  that 
C(p*(er),a)  is  monotone  increasing  for  « sufficiently  large,  and  tends  to  infinity  with 
*,  with  but  one  exception;  in  the  case  of  "linear  growth"  (f(p>  - p),  we  find  that 
C(p*(a),*)  ■ e.  However,  in  ihe  classes  of  algorithms  we  c*udy,  the  case  f(p)  «■  p 
do©i  not  arise,  provided  that  we  include  "combinatory  cost"  (see  Section  5)  in  our 
complexity  measure.  Thus  in  general,  we  find  that  the  "pointedness"  of  the  cost  curve 
near  the  minimum  increases  wfovaut  bound  as  aToo . 

Next,  we  will  show  that  for  any  f satisfying  the  hypotheses  of  Theorem  2.1,  the 
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difference  in  complexity  between  using  e method  of  fixed  order  end  a method  of 
optimal  order  tends  to  infinity  with  *. 

Proposition  2.2:  For  any  fixed  pg  < R4"*'  such  that  G'Tpg)  i.  0, 

l'metco  " C*(«)]  - +«. 

Proof:  Pick  « so  large  that  p*(c>  > pg,  and  let  pg  < p « p*(«).  If  we  write  out 
the  partial  derivative  in  the  last  term  of  (2£0),  we  *:nd  that 

dj2  C(p^)  - p”2  f(p)  e«/P  G'{p)  ♦ p-4  [e  - G(p>]  [(«  ♦ 2p)  - G(p)]  f(p)  . 

Since  pq  < p < p*(«),  we  have  G(p)  < at  it  then  follows  that  dj^p#)  I*  positive  and 
bounded  away  from  zero  as  e tends  to  infinity.  Since 

C(po^r)  - C*(«t>  - dj2  C(p,«)  [p0  - p^e)]2  / 2 
for  some  p between  p0  and  p*(e),  the  result  follows.  | 

As  for  the  cost-ratio,  a simple  calculation  stows  that 

lim.ta  C(pq^)/C*(«)  - -MX) 

in  all  of  the  examples  given  above.  Thus  there  are  a number  of  ways  in  which  we 
incur  a large  additional  cost  by  not  using  the  optimal  order. 

One  may  wonder  whether  the  result  that  optimal  order  increases  and  tends  to 
infinity  with  or  is  “reasonable."  One  way  of  determining  this  is  to  examine  actual 
numerical  tests;  we  cite  Hull  et  al  [72]  as  a well-known  example.  Since  we  are  only 
dealing  with  methods  of  fixed  order,  our  theory  decs  not  attempt  to  handle  methods 
such  as  Bulirsch-Stoer,  Krogh,  or  Gear.  However,  let  us  look  at  the  results  of  Hull  et 
al.  for  the  Runge-Kutta  methods  (which  ar£  germane  to  our  discussion—See  Section 
5.2).  Even  though  there  are  only  three  methods  (of  orders  four,  six,  and  eight)  and 
threo  error  criteria  («  - 1C-2,  10~^,  and  10”®),  Table  1 in  Huil  e!  aL  [72]  indicates  that 
the  optimal  order  does  increase  as  « decreases.  (Wa  give  mors  extensive  numerical 
data  in  Section  7.) 
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Finally,  wa  not®  that  the  restriction  that  the  grid  be  equidistant  may  b® 
w©aK®«ed  somewhat,  provided  that  we  u*e  a local  error  assure.  Indeed,  let  1 be 
partitioned  as  I «■  Ij  U-uS^,  and  now  assume  that  we  use  a grid  that  is  equidistant  on 
each  subinterval  1^,  - , I(_.  Then  the  total  complexity  is  given  by  the  sum  of  the 
complexities  of  nil  subintervals 

C\Plr-rf>L^)  8-  Z^J  CjtPj,*), 

whare  we  set 

Cj(p^i)  fj(p^)  *^p  , fj(p)  v «j(p)^p  c<p) ; 
here  Kj(p)  is  the  arror  constant  of  ?p  on  Ij.  Since  we  use  a local  error  measure,  we 
find  that  Ctp^p^e)  is  minimized  by  choosing  each  pf  to  minimize  Cj(  * , «).  Thus  the 
earlier  results  apply;  in  particular,  if  we  define  p*(a)  to  be  the  optimal  order  on  I|,  we 
fend  that  if  fj  satisfies  (2.10),  (2.11),  end  (2.15),  then  Pj*(«)  increases  and  tend2  to 
infinity  with  u. 
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Section  3 


Optimality  Within  a Basic  Sequence 


There  are  two  difficulties  with  the  approach  taken  in  Section  2.  Tho  first  has 
already  been  mentioned— we  generally  expect  the  error  coefficient  to  depend  on  the 
step-size.  The  second  is  based  on  the  fact  that  there  are  a large  number  of  p*h-order 
methods  of  a given  type,  and  we  wish  to  use  the  best  method  possible.  In  theory,  this 
would  involve  finding  a p*h-order  method  with  minimal  cost  per  step.  In  practice,  this 
is  not  often  possible*,  there  is  a gap  between  the  minimal  cost  theoretically  possible 
and  the  cost  of  the  best  method  known  So  we  now  consider  the  extension  of  the 
results  in  Section  2 to  a more  general  setting,  which  will  take  these  two  difficulties 
into  account. 

We  first  refine  our  notion  of  order.  Let  r.  ♦*!  -»  R*  be  an  error  measure, 
where  # - {*>p:  p ( Z is  a class  of  one-step  methods,  and  suppose  that  a function 
r:  R+Xl  -*  R + and  analytic  functions  , *jj  : R*  -*  R*  exist  such  that  lim  p_*g  *l(p)^p 
and  lim  xy(p)^-u  exist  and  are  nonzero  and 

(3.1)  e(*p,h)  - *(p,h)  hp  for  h « I and  p « , 

where 

(3.2)  0 < «j_(p)  S *<p,h)  S *y(p)  < -kd  for  h < 1 . 

Then  ?p  is  said  to  have  order  p with  respect  to  e,  and  $ is  said  to  be  a oasic 
sequence  (as  in  Traub  [64]h  *(p,h)  is  said  to  be  the  error  coefficient  of  *p.  (Here  we 
i .roduce  the  convention  of  attaching  the  subscripts  “l"  and  IT  to  quantities  that 
refer  to  lower  and  upper  bounds  on  complexity,  respectively.) 
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This  definition  of  order  is  similar  to  that  in  Cooper  [69]  and  Cooper  and  Verner 
72],  accept  that  we  include  a lower  bound  «|_(p)  on  n{p.h);  this  lower  bound  is 
necessary  and  sufficient  to  guarantee  that  the  order  of  a method  is  well-defined.  Note 
that  this  definition  makes  sense  for  all  values  of  h t I;  thus,  it  is  non-asymptotte  in  that 
we  do  not  require  h 1 0 in  order  for  it  to  make  sense.  Clearly,  a strong  basic 
sequence  is  a basic  sequence;  hence,  the  definition  of  order  is  an  extension  of  t*,e 
definition  of  strong  order  given  in  Section  2.  Finally,  note  that  the  order  depends  on 
the  choice  of  the  error  measure  r,  for  instance,  the  order  with  respect  to  the  local 
error  per  step  exceeds  that  with  respect  to  the  local  error  per  unit  step  by  cne. 

We  next  discuss  the  notion  of  cost  per  step.  As  pointed  out  above,  we  will 
generally  have  only  bounds  on  the  cost  c(p)  required  per  step  of  a given  pth-order 
method: 

(3.3)  Cl<p)  S c(p)  S cy(p) . 

That  is,  Cj_(p)  is  a lower  bound  on  the  minimum  possible  cost  per  step,  usually  derived 
via  theoretical  considerations,  and  Cy(p)  is  an  upper  bound  on  the  minimum  possible 
cost  per  step,  which  is  derived  by  exhibiting  an  algorithm  for  computing  ?p.  (In  what 
follows,  we  shall  assume  that  c^_ , Cy  : R+  -*  R*  are  analytic  functions.) 

We  now  wish  to  give  bounds  on  C(p,a),  the  complexity  of  finding  an  approximate 
soiution  of  (2.1)  using  the  method  such  that  ff(?p,h)  £ o-®.  Suppose  that  (2.3) 
holds.  Then  by  (3.1)  and  (3.2),  we  must  have 

(3.4)  *L(p)hp  S e~°,  i.e.,  h S hL(p,a)  *L(p)_1/pe"a/p  . 

Hence,  the  number  of  steps  n - h-*  must  satisfy 

(3.5)  n £ *l(p)^p  e°^P  ' 

Defining  (as  in  Section  2) 

<3.u) 


C(p,a)  n c(p) 
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(i.e.,  total  complexity  equals  number  of  steps  required  multiplied  by  cost  required  per 
step),  (3.3)  and  (3.5)  imply  that 

(3.7)  C(p,a)  > C^p/r)  :■  fj_<p)  e“/p  , 

where 

(3.8)  fy(p)  «l(p)^p  Cy(p). 


That  is,  regardless  of  the  algorithm  used  to  compute  ^p,  the  total  complexity  of  finding 


an  approximate  solution  of  (2. 1 ) must  exceed  C^p,*). 

On  the  other  hand,  we  find  that  in  order  to  use  ?p  to  find  such  an  approximate 
solution,  it  suffices  (by  (3.1)  and  (3.2))  to  take 

(3.9)  Ky{p)  hp  - e~a,  i.e.,  h « hy(p,«)  :■»  ityfpr^P  e"a/p  . 
so  that  we  need  only  take  n steps,  where 

(3.10)  n « «y(p)l/P  e°/P  . 

(As  in  Section  2,  the  value  of  n given  by  (3.10)  need  not  be  an  integer;  again,  this  is 
handled  as  in  Traub  and  Wofniakowski  [76].)  Thus  (3.3),  (3.6),  end  (3.10)  imply  that 

(3.11)  C(p,a)  S Cy(p,e)  :»  fy(p)  ec/P  , 

where 

(3.12)  fy(p)  :=  xy(p)^P  Cy(p). 

That  is,  there  exists  an  algorithm  for  computing  pp  such  that  the  total  complexity  of 
finding  an  approximate  solution  of  (2.1)  equals  Cy(p,a).  We  summarize  the  above 
results  in 
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Theorem  3.1:  Let  C(p,«)  be  the  complexity  of  finding  an  approximate  solution  of 
(2.1),  using  the  method  *p,  with  *(*p,h)  £ e"®.  Then 

(3.13)  C^(p,or)  £ C{p,o)  £ C(j{p,o)  | 

where  Cy  and  Cy  are  given  by  (3.7)  and  (3.11).  Moreover,  if  h - h(p^r)  is  the  maximal 
step-size  for  the  method  ?p  such  that  *(?p,h)  £ e-®,  then 

(3.14)  hy(p,or)  £ h(p,«)  £ hy(p,or) . | 

Next,  we  consider  the  problem  of  optimality.  Define  the  optimal  complexity  by 

(3.15)  C*(«)  :«  inf  {C(p,er>:  ?p  € *} . 

We  are  interested  in  bounds  for  C*(a).  These  are  derived  in 

Lemma  3J.:  Let  fy  and  fy  satisfy  (2.10)  and  (2.11),  and  suppose  that  fy  and  fy 
respectively  yield  (via  (2.13))  Gy  and  Gy  satisfying  (2.15).  Then  Gy  and  Gy  have 
respective  inverse  functions  pL*,  py*:  R++  -»  R++  such  that  for  all  p ( R*4, 

(3.16)  CL*(«)  CL(pL*<«),«)  £ CL(p,«) 

and 


(3.17)  Cy*(«)  Cy(py*(«),«)  £ Cy(p,a)  , 

with  equality  in  (3.16)  (respectively,  (3.17))  if  and  only  if  p - py*(a)  (respectively,  p - 

Proof:  This  is  an  immediate  corollary  of  Theorem  2.1.  | 

We  call  py*(a)  (respectively,  py*(o))  the  lower  (upper)  optimal  order.  Cy*(a) 
(respectively,  Cy*(a))  the  lower  (upper)  optimal  complexity,  and 

(3.18)  hy*(«)  hj  (py*(a),o)  (respectively,  hy*(a)  hy(py*(a),a)) 

the  lower  (upper)  optimal  step-size.  Combining  (3.13),  (3.15),  and  Lemma  3.1,  we  have 
Theorem  3_.2:  Let  fy  and  fy  be  as  in  Theorem  3.1.  Then 
CL*<«)  £ C*(«>  £ Cy*(«)  . | 

Note  that  if  we  define  p*(a)  by 
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C(p*(cr),a)  - C*(«) , 

wo  can  make  no  statement  relating  p*(«),  py*(«),  and  py*(«).  This  is  because  we  only 
have  bounds  for  C(p,a)j  we  do  not  know  C(p,a)  itself.  In  fact,  it  is  important  to  realize 
\ 'hat  py*(«)  and  py*(«)  tell  us.  First,  consider  py*(«).  We  can  achieve  a complexity  of 
Cy*(«)  by  using  a step-size  of  hy*<«),  along  with  the  method  of  order  py*(«).  This  will 
give  optimal  complexity  within  the  sequence  of  algorithms  for  computing  4,  with  cost 
per  step  of  given  by  cy{p).  Naxt,  consider  py*(a).  It  is  of  perhaps  theoretical 
rather  than  computational  interest,  in  that  we  cannot  compute  with  it.  What  does 
interest  us  is  Cy*{<r),  since  it  limits  the  theoretical  improvement  in  Cy*{«).  Thus,  we 
are  interested  in  py*(«)  solely  as  a means  of  computing  Cy*(«). 

We  now  consider  behavior  of  these  quantities  as  a increases  and  tends  to 
infinity. 

Theorem  3.3:  Let  fy  and  fy  be  as  in  Theorem  3.1.  Then 

(1.)  py*(a),  Py*{«),  Cy*(«),  and  Cy*(*>  increase  monotonically  and  tend  to 

infinity  with  «. 

(2.)  If  there  exists  an  My  > 0 such  that  *y(p)^p  £ My  for  all  p,  then 
lim  infgfog  hy*(«)  > 0 if  «/py*{«)  is  bounded  as  atoo. 

(3.)  If  there  exists  an  My  > 0 such  that  jty(p)l/p  * My  for  all  p,  then 
lim  inf hy*(a)  > 0 only  if  a/py*(a)  is  bounded  as  otoa. 

Proof:  To  prove  (1.),  it  suffices  to  apply  (1.)  and  (2.)  of  Theorem  2.2  to  py*  and 
Cy*,  and  to  py*  and  Cy*.  The  proof  of  (2.)  and  (3.)  is  similar  to  the  proof  of  (3.)  in 
Theorem  2.2.  | 

Note  that  (1.)  in  Theorem  3.3  does  noj,  state  how  p*(a)  varies  with  a;  as  we  have 
pointed  out  above,  no  statement  about  p*(a)  may  be  obtained  from  the  information 
available.  However,  it  is  easy  to  see  that  C*(a)  increases  monotonically  with  a and  that 

limotooc*<«)  “ +0°  • 
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Thus,  we  have  extended  the  optimality  theory  of  Section  2 to  a more  realistic 
situation.  In  Sections  5 and  6,  the  techniques  of  this  section  will  be  applied  to  some 
important  basic  sequences  of  one-step  methods;  we  will  see  that  the  conclusions  of 
Lemma  3.1  and  Theorems  3.2  and  3.3  hold  for  these  basic  sequences. 


wm, 


Section  4 

Normality  and  Order-Convergence 


Let  4 be  a basic  sequence  with  respect  to  the  error  measure  r,  we  say  that  4 is 
order-convergent  if  there  exists  an  tig  > 0 such  that 

(4.1)  timptoo  *l/P>  hP  “ 0 for  h£  hg  . 

Clearly,  the  order  convergence  of  ♦ implies  that  limpf^  e(*p,h)  - 0 for  h £ hg.  We 
use  the  term  "order-convergence"  rather  than  "convergence,"  since  the  latter  term 
appears  extensively  in  the  literature  (e.g.,  Henrici  [62])  and  is  always  used  to  mean  a 
"step-size  convergence,"  i.e.,  lim^  r(fth)  « 0 for  a fixed  method  ?. 

It  is  intuitively  plausible  that  as  the  order  of  an  approximation  increases,  the 
approximation  should  improve,  especially  when  one  is  trying  to  approximate  a very 
smooth  function.  Unfortunately,  Gear  [71]  points  out  that  an  increase  in  order  need 
not  always  decrease  the  error.  This  situation  appears  in  other  situations  in  numerical 
mathematics;  for  instance,  the  family  of  Newton-Cotes  quadrature  formulae  is  no]: 
order -convergent.  But  suppose  there  exists  a step-size  hg  > 0 for  which  the  upper- 
bound  error  is  exponentially  bounded  for  p sufficiently  large;  that  is,  there  exists 
A > 0 and  pg  € 1*+  such  that 

(4.2)  k^Kp)  hgP  £ AP  for  p > Pg  . 


If  we  define 

My  :«  max  { max1SpSp  fri/p)1^} , Ahg-1} , 

we  then  have 
(4.3) 


*(*p,h)  £ (Myh)P  for  h £ hg  , p € Z . 
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Note  that  the  bound  in  (4.3)  is  similar  to  that  given  by  Cauchy’s  Integral  Theorem 
(Ahlfors  [661  pg.  122)  on  the  normalized  derivatives  of  an  analytic  function.  In  fact, 
for  several  classes  of  methods,  the  bound  (4.3)  holds  whenever  the  solution  of  (2.1)  is 
analytic. 

We  also  formalize  a weakened  version  or  (4.3),  which  will  be  important  in  our 
study  of  one-step  methods.  Let  * be  a basic  sequence,  and  suppose  that  for  each 
(xq,v)  € SW,  there  is  a sequence  {hp:  p € Z ++}  c I and  a positive  constant  My  such 


that 


(4.4)  #(*p,h)  <;  (Myh)P  if  h s hp  ; 

then  # is  said  to  be  normal.  Note  that  (4.3)  implies  (4.4),  while  (4.4)  implies  (4.3)  only 
when  the  sequence  {hp}  has  non-vanishing  support: 

(4.5)  hg  :■*  lim  intp^  hp  > 0 . 

If  h^  » 0,  normality  gives  an  exponential  upper  bound  on  the  sequence  of  principal 
error  functions  (Section  3.3-5  of  Henrici  [62]),  which  are  an  asymptotic  measure  of  the 
error  as  h A 0. 

There  is  a simple  relation  between  normality  and  order-convergence. 

Proposition  4.1:  * is  order-convergent  if  and  only  if  * is  normal  with 
nonvanishing  support. 

Proof:  If  (4.1)  holds,  then  (in  particular)  we  have  lim^o,  *y(p)  hgP  - 0,  so  that 
«y(p)  hQP  S 1 for  p sufficiently  large;  i.e.,  (4.2)  holds  with  A ■ 1.  Then  (as  in  tho 
discussion  above)  (4.3)  holds,  implying  normality  with  finite  support. 

Conversely,  if  (4.4)  holds  with  finite  support,  we  pick  a positive  h^  which  is  less 

than 

f min  {My-*,  inf  {hp:  p ( Z++}  } . 


(Not;  that  n > 0 by  (45).)  Let  h s hQ  be  given,  so  that  tor  some  I with  0 < I < 1,  we 
have  h - (1  - if  we  define  ay  by 

«y(p)  My°, 

we  find  (since  h < hp)  that 

e(^p,h)  d (Myh)P  . 

Thus 

*y(p)  hP  - (Myh)P  - (My  (1  - !)*)*  S (1  - J)P 
(the  last  step  since  f £ My"*),  so  that  (4.1)  holds,  d 

We  are  now  interested  in  normality  and  order-convergence  for  a specific  error 
measure  r,  we  will  be  interested  in  cyy , *j_ , and  cq  , which  are  (respectively)  defined 
to  be  the  maximum  local  error  per  unit  step,  local  error  per  step,  and  global  error  per 
step  over  the  grid.  It  is  easy  to  see  that  a normal  (order-convergent)  sequence  * - 
{*p;  p « Z**}  with  respect  to  naturally  yields  a normal  (order-convergent) 
sequence  * - p;  p € Z ++}  with  respect  to  ayy  by  setting  f p *p+1  for  p ^ Z 

We  now  look  at  the  relationships  between  ryj  and  *g. 

Proposition  4.2;  Let  v have  Lipschitz  constant  K on  F^,  and  let  ♦ be  normal 
(respectively,  order-convergent)  with  respect  to  cyj,  with  My  in  (4.4)  Independent  of 
Xq  « domain(v).  Then  * is  normal  (respectively,  order-convergent)  with  respect  to  #g. 

Proof;  Let  p be  the  exact  relative  increment  function  of  (2.1)  (as  defined  in 
Henrici  [62]),  so  ‘hat 

x(ti+i>  - x(tj)  ♦ h p(x(tj),  h ) . 

Subtract  (2.2)  (with  p replaced  by  /p)  from  the  above  to  get 

»i+l  " «i  ♦ h [p(x(tj),h)  - ^p(xj,h)] , 
where  ej x(tj)  - X|  for  C £ i £ n.  Thus 
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n®i+lli  s ||»ill  ♦ h ||/Kx(tj),h)  - p(Xj,h)||  * h ||p{xith)  - ^p<X|,h)!| 

S (1  ♦ hK)  !|sj||  + MyP  hP+1  if  h £ hp  ; 

this  last  step  follows  from  the  Upschitz  condition  and  the  "uniforsi"  normality  with 
respect  to  #li>  By  Lemma  1.2  of  Henrici  [62]  and  the  condition  Sq  « 0,  we  have 
||e;||  £ IT1  fU  + hK)1  - 1]  (MyhjP 
£ IT1  ((I  * hK)n  - 1]  (MyhJP 
S K"1  <eK  « i)  (MyhjP 

for  ail  ij  this  gives 

.fQtep.h)  £ IT1  (eK  - t)  (MyhjP  £ (Mh)P  if  h £ hp> 
for  a suitably-defined  M > 0.  This  proves  the  normality  part;  the  remainder  of  the 
result  follows  from  Proposition  4.1.  | 

If  it  is  undesirable  to  use  the  "uniform  normality"  (i.e,  the  condition  that  My  be 
independent  of  xq  i domain(v)  in  (4.4)),  we  may  use  the  following  result. 

Proposition  4.3:  Let  v be  Upschitz  continuous,  let  * be  normal  (respectively, 
order-convergent)  with  respect  to  ^0  a*»d  suppose  that  there  exists  a X > 0 such  thst 
for  all  ipp  < $ and  all  x,  y « R^, 

l!*p<x)-*p(y)H  s Xpl|x-y||  . 

Then  w is  normal  (respectively,  order-convergent)  with  respect  to  *q. 

Proof:  Immediate  from  Theorem  3.3  of  Kenrici  [62] . | 

Thus  normality  for  *q  follows  from  normality  for  oy  , a Upschitz  condition  on  v 
and  the  elements  of  #,  and  a linear  upper  bound  on  the  Upschitz  constants  for  the 
elements  of  ♦. 

We  now  discuss  the  problem  of  finding  uniform  lower  bounds  on  the  error  which 
are  similar  to  the  uniform  upper  bounds  which  normality  provides.  This  will  amount  to 
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restriction  of  the  admissible  problem  class  so  as  to  guarantee  thot  the 
problems  are  "sufficiently  difficult.”  However,  this  restriction  may  be  abandoned  if  we 
are  interested  only  in  upper  bounds.  Wg,  shall  assume  throughout  the  rest  o£  this 
paper  that  there  ii  an  M-^  > 0 (which  wiil  generally  dtsperd  on  #,  r,  and  the  problem 
(xq,v))  such  that 

. (4.d}  c(*p,h)  ;>  (M^  for  h < I . 

Note  that  (4.6)  will  hold  for  nrsy  situation  in  which  th..  re  is  no  or  dor -convergence,  or  in 
which  the  order-convergence  (if  any)  is  no  faster  than  an  exponential  decay; 
moreover,  in  the  methods  we  consider  in  Sections  5 and  6,  (4.6)  is  a consequence  of 
the  assumption  that  all  derivatives  assume  the  (sharp)  worst-case  upper  bound 
provided  by  Cauchy’s  estimate.  It  is  Hear  that  if  (4.6)  holds  for  , it  holds  for  *l«j  » 
if  (4.6)  holds  for  r^j  and  if  the  matrix  Vp  has  only  non-negative  entries  (with  at  least 
one  positive  entry),  then  (4.6)  holds  for  #q  . 

It  is  possible  to  present  a simplified  version  of  the  expressions  derived  in 
Section  3,  under  the  assumption  that  ♦ is  order-convergent.  We  first  look  at  the 
complexity  of  a single  method  within  an  order -convergent  basic  sequen  e. 

Theorem  4.1:  Let  * be  order-convergent  with  respect  to  *.  Then 
CL(p,*)  S C(p^r)  S Cy(p,cr) , 

where 

CL(p^t)  :-  Mj_  c^(p)  ea/P  and  C^p^t)  :-  My  CyCp)  e*/p  . 

Proof:  This  is  an  immediate  corollary  of  Theorem  3.1  and  the  definition  of 
order-convergence.  1 

We  may  now  do  the  optimality  theory  of  Section  3,  finding  that 

(4.7)  G^(p)  - p2  cl^pVc^p)  and  Gj/p)  - p2  cy'^/c^p)  . 


Note  that  the  assumptions  (2.10)  and  (2.11)  now  state  that  c^(p)  and  Cy(p)  must  be 
positive  for  p > 0 and  tend  to  infinity  with  p,  which  is  a natural  way  to  expect  the  cost 
per  step  to  behave.  The  results  stated  in  Theorems  3.2  and  3.3  hold  as  before. 
Moreover,  it  should  be  noted  that  the  My  and  needed  in  (2.)  and  (3.)  in  the 
statement  of  Theorem  3.3  are  precisely  the  My  and  in  (4.4)  and  (4.6).  Thus 
lim  inf^oj  hy*(a)  > 0 if  o/py*(«)  is  bounded  as  « t oo,  «ivd  afp^(a)  is  bounded  as 
a T oo  if  iim  inf ^3,  hL*(a)  > 0. 

Thus,  the  order-convergence  of  a basic  sequence  is  useful  in  simplifying  the 
analysis  of  its  complexity.  Of  the  three  basic  sequences  we  will  study  in  this  paper, 
two  are  Known  to  be  order-convergent.  The  proof  of  the  order-convergence  of  the 
class  of  Taylor  series  methods  is  a simple  consequence  of  the  Cauchy  estimate;  that  of 
the  order-convergence  of  the  (non-op^aUy  ordered)  nonlinear  Brent-Runge-Kutta 
methods  (given  in  Appendix  B)  invcives  using  some  classical  results  u.i  orthogonal 
polynomials  to  sharpen  the  proof  a Brent  [74]  . (We  note  that  it  is  not  Known 

whether  the  cptimally-ordersd  ->:„‘;..ear  Brent-Runge-Kutta  methods  ere  order- 
convergent;  it  does  appear  liKeiy  that  they  are  normal  with  vanishing  support. 
However,  we  do  not  pursue  this  class  of  methods,  because  of  their  high  combinatory 
cost,  as  indicated  in  Section  6.) 

It  is  not  Known  whether  the  linear  Runge-Kutta  methods  found  in  Cooper  [69] 
and  in  Cooper  and  Verner  [72]  are  order-convergent;  the  best  result  Known  is  the 
(My  logtp+e))13  result  given  in  Appendix  A,  which  involves  strengthening  the  original 
proof  with  other  estimates  from  the  theory  of  orthogonal  polynomials.  But  it  should  be 
pointed  out  that  there  does  exist  a class  of  order-convergent  linear  Runge-Kutta 
methods;  this  is  the  sequence  given  by  using  the  weights  .and  abscissae  for  Gauss 


29 


SisBlfeJSS 


quadrature  in  the  methods  defined  on  page  144  of  Stetter  [73J  The  problem  with  this 
class  of  methods  is  that  each  step  of  requires  2_p(pvl)l  function  evaluations;  the 
prohibitive  cost  per  step  outweighs  by  far  any  advantage  to  be  gainod  from  the 
order-convergence.  Thus,  the  question  of  whether  there  exist  any  order  convergent 
linear  Runge-Kutta  methods  which  are  more  efficient  (he.,  have  smaller  cost  per  step) 
remains  open. 


Section  5 


mmz 


Applications  to  Systems  of  Differential  Equations 


In  the  next  two  sections,  we  apply  the  theory  developed  in  the  proceeding 
sections  to  two  of  the  most  commonly-used  classes  of  methods,  l.e^  Taylor  series 
methods  and  Runge-Kutta  methods.  In  Section  5,  we  shall  treat  the  complexity  of 
systems  of  differential  equations,  i.e.,  problems  of  the  form  (2.1)  for  which  v is  an 
operator  on  R^,  wnere  N is  an  arbitrary  positive  integer.  In  Section  6,  we  shall 
restrict  cur  attention  to  the  scalar  case,  i.e^  the  case  where  consists  of  functions 
v:  R -*  R;  for  this  case,  Brent  [74]  has  discovered  a class  of  "nonlinear  Runge-Kutta 
methods." 

Before  discussing  the  complexity  of  these  basic  sequences,  we  fix  our  error  and 
cost  measures.  For  the  sake  of  definiteness,  we  shall  choose  »q  as  our  error  measure; 
that  is,  we  willl  be  interested  in  the  global  nrror,  rather  than  the  local  error  per  step 
or  per  unit  step.  However,  the  other  error  measures  may  be  used  with  a slight 
modification  of  the  discussion  contained  in  the  sequel 

We  now  make  precise  cur  notion  of  cost  We  win  bg  concerned  with  the  total 
number  cf  arithmetic  operations  required.  Let  * be  a given  basic  sequence.  As  in 
Traub  and  Woiniakowski  [76],  we  shall  express  the  cost  per  step  ^sociated  with  in 
the  form 

(5.0.1)  c(p/ :«  e(J2p(v))  + d{p)  . 

Here  S£p(v)  <s  the  information  about  v required  to  perform  one  step  of  *p,  and  we 
write  e(9Jp(v))  for  the  informational  cost  of  *p;  we  call  (ftp)  the  combinatory  cost 
Of  |?p- 


Note  that  we  explicitly  indicate  the  dependence  of  y)p  on  v,  so  that  we  may 

compare  the  cost  of  (say)  an  evaluation  of  v with  a scalar  arithmetic  operation. 

Basically,  e(52p(v))  measures  the  cost  of  getting  new  data  about  v required  by 

while  d(p)  measures  the  cost  of  combining  this  new  data  to  get  an  approximate  value 

of  the  solution  at  a new  point.  For  example,  Euler's  method  in  R^ 

xi+l  “ *i  + Mxp 
N 

has  informational  cost  e(vj) , where  vj  , ~ , v^  are  the  components  of  v and  for 
any  function  u:  R^  -*  R,  we  define 


(5.0.2) 


e(«)  :•  cost  of  evaluating  <*  at  one  point  . 


The  combinatory  cost  is  2N  arithmetic  operations,  he*  one  scalar  multiplication  and  one 
scalar  addition  for  each  of  the  N components. 

Finally,  we  now  assume  that  0 and  t?  have  been  chosen  so  as  to  guarantee  that 
the  solution  of  (2.1)  is  analytic  on  l Thus  Cauchy's  Integral  Theorem  guarantees  the 
existence  of  a positive  M such  that  for  all  positive  integers  p,  we  have 

(1/p!)  H*(p)(t)||  S MP  for  t U . 
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5.1.  Taylor  Sari**  Methods 

The  class  #-j*  of  Taylor  series  mettods  is  defined  by  expanding  x in  a truncated 
Taylor  series.  Thus  the  increment  function  ?p  is  given  by 

(5.1.1)  *p(xj,h}  2^  v(K)(Xj)  hk  / (K+l)! , 

where 

(5.1.2)  v(W(xj)  <d/dt)k  [v(x(t))]  | x(t)  . Xj  . 

The  usual  method  of  computing  (5.1.2),  as  described  in  "classicar  numerical  analysis 
texts  such  as  Henrici  [621  invokes  the  chain  rule.  This  quickly  leads  to  expressions  of 
horrifying  complexity;  for  this  reason,  most  texts  quickly  abandon  the  discussion  of 
high-order  Taylor  series  methods. 

We  are  interested  in  faster  algorithms  for  computing  *p.  First,  we  address  the 
problem  of  a lower  bound  for  the  combinatory  cost  d(p). 

Proposition  5.1.1:  There  exists  a constant  > 0 such  that  any  sequence  of 
algorithms  for  computing  $ y must  satisfy 

(5.1.3)  d(p)  S aLpN  . 

Proof:  Any  algorithm  for  computing  $>p  requires  the  information 
S2p(v)  {D^v:  0 £ |/Jj  S p - 1}  . 

(We  use  the  standard  multi-index  notation  found  in  Friedman  [69P  It  is  then  easy  to 
see  that  the  above  set  has  C\pty  (as  p T as)  distinct  elements,  which  are  (generally) 
independent;  this  is  an  immediate  consequence  of  Problem  1 1 in  Chapter  I of  P6lya  and 
SzegB  [251  Thus  (5.1.3)  gives  a linear  lower  bound,  g 

ftote  that  the  constant  in  (5.1.3)  depends  on  N.  Since  we  are  treating  the 
case  where  N is  fixed  and  p is  allowed  to  very,  wa  will  not  indicate  this  dependence 
explicit!/.  We  now  sea  how  close  we  can  get  to  on  optimum  value  for  d(p). 
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Theorem  5.1.1:  There  exists  a constant  ay  > 0 such  that  the  combinatory  cost 
d(p/  of  computing  ^p  € satisfies  the  bound 

(5.1.4)  d(p)  £ ay  p^  In  (p+e)  . 

Proof:  We  first  consider  the  case  N ~ 1.  Note  that  x(h)  is  the  zero  of 

(5.1.5)  F(z)  :«  SI  d*  / v<*)  - h • 

*0 

As  in  Brent  and  Kung  [76],  we  consider  the  formal  power  series 

P(s)  F(xq+s)  - F(x0) , 

where  s is  an  indeterminate.  Let  V be  the  power  series  reversion  of  P.  Adopting  the 
notation  of  Brent  and  Kung  [76],  we  see  that 

x(s)  - X0  + V(s)  - X0  + Vp(s)  + 0(sP+1)  . 

By  the  uniqueness  of  the  Taylor  coefficients  of  an  analytic  function,  we  see  that 

*>p(x<)>h)  “ h_1Vp(h) . 

Since  the  number  Vp(h)  can  be  computed  in  0{p  In  p)  operations  from  the  Taylor 
coefficients  of  v (by  Theorem  6.2  of  Brent  and  Kung  [76]),  the  result  for  N - 1 follows. 

For  N > 2,  we  use  Newton’s  method  (Rail  [69])  applied  to  the  formal  power 
series  operator  P given  by 

(PyHs)  :*  y(s)  - x0  - Jq  v(y(r))  dr  ; 

clearly,  the  formal  power  series  x(s)  is  the  zero  of  P.  The  algorithm  itself  is  defined 
recursively.  Let  a formal  power  series  X(p)(s)  satisfying 

x(p)(s)  - x(s)  + 0(sp+1) 


be  given.  Precompute 

(5.1.6)  w(s)  :«  Jq  v(X(p)(r))  dr  - x0  - X(p)(s)  + 0{s2p+2) , 

(5.1.7)  Q(s)  :-  W(x(p)(s))  + 0(s2p+2) , 
and  let  U(q)(s)  :«  0.  Then  set 

x(2p+l)<s>  *•-  x(p)(s)  + U(P+1)(S)  ’ 


m 
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where 

(5.1.8)  U(k+|)(s)  :-  j*Q  Q(f)  u^)M  dr  + w($)  + (Xs2P+2) , 0 S K 5 p . 

Following  the  proof  giver,  in  Rail  [691  we  find  that 

x(2p+i)(s)  - x(s)  + 0{s2P+2)  . 

We  need  only  consider  the  cost  T(p,N)  of  computing  the  series  *(p)(s)  in 
determining  d(p),  since  x(h)  may  be  recovered  from  the  formal  power  series  in  CXp) 
operations.  Clearly,  we  have  the  recursion 

(5.1.9)  T(2p+1,N)  S T(p,N)  + Tg  + T;  + Tg  , 

where  Tm  is  the  cost  of  step  (5.1.m)  for  m - 6,  7,  8.  Let  COMP(p,N)  be  the  time 
required  to  find  the  first  p terms  of  the  formal  power  series  f(y^(s),  ... , y^fs)),  where 
f,  Yi.  ~ i yfyj  are  formal  power  series,  and  y^,  ...  , y»g  have  zero  constant  term. 
Theorem  7.1  of  Brent  and  Kung  [76]  states  that 

C0MP(p,2)  - 0(p2  In  p) , 
and  it  is  easy  to  show  that  for  any  N * Z ++J 

COMP(p,N+i)  - (Xp  COMP(p.N))  . 

Thus  for  N £ 2,  we  have 

(5.1.10)  C0MP(p,N)  - 0{pNlnp), 
and  so  we  see  that 

T6  + T7  - CK(2p+l)N  In  p) . 

Finally,  let  MULT(p)  be  as  in  Brent  and  Kung  [761  we  see  that 

Tg  - (p+1)  [N2  MULT(2p+l)  + 0(p)]  - (X(2p+1)2  In  p) 
if  Fast  Fourier  Transform  multiplication  (Borodin  anj  Munro  [75])  is  used.  Since  N £ 2, 
we  have 

(5.1.i  d T6  + T?  + Tg  . 0«2p+l)N  In  p) , 

and  so  (5.1.9)  and  (5.1.11)  imply  that 


r ■>  1 -■  (v 


T».j: 
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T<p,N)  - CKpN  In  p) , 

which  completes  the  proof.  | 

(Note  that  the  second  algorithm  is  inferior  to  the  first  algorithm  when  applied  to 
the  scalar  case  N - 1,  where  we  find  that  the  second  algorithm  requires  CKp^  In  p) 
arithmetic  operations.) 

We  now  determine  bounds  on  C(p,a).  First,  consider  lower  bounds.  Clearly, 
there  exists  ej_(v)  2 0 such  that 

(5.1.12)  e(D^Vj)  2 eL(v)  (lsi£n,|0|€Z+)  . 

Since  $lp(v)  has  O(pN)  elements,  there  exists  a constant  b^^  > 0 such  that 

(5.1.13)  e(fflp(v))  2 bL  eL(v)  pN  . 

From  (5.1.3)  and  (5.1.13),  we  have  a lower-bound  cost  per  step  of 

(5.1.14)  cL(p)  - |>l  + bL  ei_(v)]  pN  . 


This  leads  to 

Theorem  5.1.2:  C^p^r)  - taL  + eL^v^  PN  e°^P  • 

Proof:  This  is  an  immediate  consequence  of  (4.6)  and  (5.1.14).  | 

Note  that  f|_(p)  :-  M^c^p)  satisfies  the  conditions  of  Theorem  3.2.  So  the 
optimality  theory  of  Section  3 holds.  In  particular,  we  have 

Theorem  5,1.3:  Cl*(o)  - [a^  + b^  e^v)]  (e/N)^  oN  . 

Proof:  From  (4.7)  ^ and  (5.1.14),  we  find  that  G[_(p)  - Np,  so  that 
p [_*(«)  - q/N  and  h^to)  - (M^e^)"!  . 

The  result  follows  by  letting  p - Pl*(«)  in  the  definition  of  C^a).  1 

However,  recall  that  we  assumed  that  the  non-identical  mixed  partial  derivatives 
of  v are  independent.  There  are  a number  of  systems  for  which  this  is  not  true  (for 
instance,  constant  coefficient  linear  systems);  for  such  systems,  it  is  clear  that  we  may 


36 


be  able  to  use  the  extra  information  of  non-independence  to  find  algorithms  that  are 
faster  than  the  lower  bounds  given  above.  However,  we  will  ignore  this  case  and  only 
consider  the  problem  for  a "general"  function  v. 

Next,  we  turn  to  upper  bounds  on  the  complexity.  Theorem  5.1.1  tells  us  how  to 
combine  the  necessary  information  to  get  the  solution  at  a new  grid-point;  we  need 
only  measure  the  cost  of  getting  the  information.  So,  let 

e(K)(v)  - max  {e{D*Vj):  1 2 i £ ty,  M - k}  . 

Using  the  result  in  P6iya  and  SzegS  [25],  we  see  that 

(5.1.15)  e(flp(v))  2 N e(k)(v>  (N+k-l)l  / [kl(N-l)!]  . 

Unfortunately,  the  right-hand  side  of  (5.1.15)  does  not  fit  our  general  model,  so  we 
must  assume  that  we  know  hew  e^(v)  changes  as  k increases.  We  win  consider  the 
case  where  the  cost  of  derivative  evaluation  is  bounded:  that  is.  we  will  assume  that 

(5.1.16)  e^(v)  S ey(v) 

for  some  ey(v)  independent  of  k.  Other  cases  (e.g.,  e^(v)  « 0(km)  for  some  m > 0) 
may  be  analyzed  in  a similar  manner;  of  course,  they  will  give  different  results.  By 
(5.1.15)  and  (5.1.16),  there  is  a by  > 0 such  that 

(5.1.17)  e(g?p(v))  S by  e(j(v)pN . 

From  (5.1.4)  and  (5.1.17),  we  have  an  upper-bound  cost  per  step  of 

(5.1.18)  Cy(p)  - ay  pN  In  (p+e)  + by  ey(v)pN  . 

This  leads  to 

Theorem  5.1.4:  There  exists  an  My  > 0 such  that 

Cy(p/»)  - My  [ay  pN  In  (p+e)  + by  ey(v)pN]  e“/p  . 

Proof:  By  Cauchy’s  Integral  Theorem,  there  exists  a B > 0 such  that 

|||x(k+1)|||  / (k+1)!  5 Bk  , 

where  we  define 
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(5.1.13)  lllylll  :-  max  t € l ||y(t)|| 

for  any  y:  I -»  R^.  Thus  by  Section  3.3-3  of  Henrici  [«21  we  see  that  a Lipschite 
constant  for  *p  in  *y  is  given  by 

Z*-0  |Hx(k+1>IH  hK  / (k+l)!  5 SK^0(Bh)K  5 Lt-U-Bhor1, 
provided  that  h S h^  < B"*.  By  Section  3.3-4  of  Henrici  [62]  and  Proposition  4.3,  there 
exists  an  My  > 0 such  that 


*G(*p,h)  S (Myh)P. 

The  resuit  now  follows  from  Theorem  4.1  and  (5.1.18).  | 

We  are  now  ready  to  consider  the  optimal  p for  Cy(p,«). 

Theorem  5.1.5: 

(1.)  For  all  a > 0,  there  exists  py*(«)  such  that  (3.17)  holds. 
(2.)  Py*(«)  increases  monotonically  with  a,  and 

Py*(«)  ~ «/N  as  « t oo . 

(3.)  Cy*(«)  increases  monotonically  with  a,  and 

Cy*(«)  ~ My  ay  (e/N)N  «N  In  « as  « t co . 

(4.)  hy*(«)  ~ (My  eN)_1  as  « tea. 

Proof:  Clearly  cy  satisfies  (2.10)  end  (2.11).  Now  write 

6y(p)  ■ Gj(p)  ♦ G2(p) , 

where 

G^p)  - N p and  G2(p)  - * p2/D2<p) ; 

here  we  set 


□2(p)  :«  (p+e)  [<p+«)  In  (p+e)  ♦ 1]  and  p ay  / [by  ey(v)]  . 

We  see  immediately  that  G}  satisfies  (2.15)?  a straightforward  calculation  shows  that 
G2'(p)  - p [D(p>]"2  (rp  [In  (p+e)]  - i]  + 2e[r  In  (p+e)  + 1]} , 
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so  that  G2/<p)  > 0 for  p > 0.  Thus  G£  satisfies  (2.15),  which  shows  that  Gjj  satisfies 
(2.15).  Hence  py*  and  Cy*  behave  as  described  in  Theorem  3.3. 

Since  py*(«)  goes  to  infinity  with  a,  we  see  that 

« - Gy(py*(«))  ~ N Py*(«)  + Py*(«)  / In  Py*(«)  ~ N Py*(«)  , 
which  gives  the  asymptotic  estimate  in  (2.).  The  rest  of  the  Theorem  follows  from  this 
estimate.  I 

Unfortunately,  the  estimates  given  above  are  only  asymptotic  as  a t oo;  this  will 
be  typical,  since  many  of  the  equations  to  be  solved  involve  products  of  logarithmic 
and  polynomial  terms,  and  thus  cannot  be  solved  exactly.  On  the  other  hand,  these 
asymptotic  expressions  are  sufficient  for  our  purposes,  since  they  describe  how 
quickly  py*(«)  and  Cy*(«)  increase  with  «. 

Note  that  as  a tends  to  infinity,  Cy*(*)  becomes  independent  of  ey(v),  which 
measures  how  hard  it  is  to  evaluate  the  derivatives  of  vj  this  is  because  the 
combinatory  cost  eventually  overwhelms  the  informational  cost.  This  kind  of  buhavior 
will  be  typical  of  the  complexity  analyses  in  this  paper.  Finally,  note  that  the  bound 

(5.1.20)  CL*<«)  - 0(«N)  S C*(«)  S 0(«N  In  e)  - Cy*(a)  as  « t co 
implies  that 

Cy*(«)  / Cj>)  * 0(ln  or)  as  a T oo  ; 

this  indicates  the  gap  in  our  knowledge  of  the  complexity  of  solving  (2.1)  via  Taylor 


series  methods. 
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S.2  Linear  Runge-KuNa  Method* 

For  many  functions  v,  caculation  of  the  derivatives  required  by  Taylor  series 
methods  is  prohibitively  expensive.  For  this  reason,  we  are  interested  in  methods 
which  use  Information  that  is  somewhat  more  readily  available.  In  particular,  we  will 
consider  methods  that  use  only  evaluations  cf  v,  combined  in  a highly  structured 
manner.  We  say  that  ^lRK  *s  a c*ass  °*  *lnaar  Runge-Kutta  methods  (abbreviated,  "IRK 
methods")  if  each  increment  function  ?p  may  be  written  in  the  form 

(5.2.1)  *p(xith)  Zj^1  X,,  K, 

where 

(5.2.2)  Kj  :*  v(Xj  + h Z^jj  Xjj  kj)  for  0 S I S s - I , 

the  integer  s - s(p)  is  said  to  be  the  number  of  stages  of  *p;  the  number  of  stages  is 
equal  to  the  number  of  times  the  vector  function  v must  be  evaluated.  (In  order  to 
simplify  notation,  we  will  not  explicitly  indicate  the  dependence  of  Xjj  and  Kj  on  p.)  The 
method  y»p  defined  by  (5.2.1)  and  (5.2.2)  is  explicit  in  that  kj  depends  only  on 
Kq  , - , k|_^5  see  Butcher  [64a]  for  a discussion  of  semi-cxplicit  and  implicit  method*. 

Since  the  function  is  (in  practice)  always  evaluated  by  using  the  obvious 
algorithm  suggested  by  its  definition,  we  shall  identify  an  algorithm  for  evaluating  *>p 
with  itself.  Thus  the  problem  of  finding  the  best  algorithm  for  evaluating  *p  in 

<&LRK  is  equivalent  to  the  problem  of  finding  the  best  basic  sequence  of  LRK  methods 
possible.  This  is  related  to  the  problem  of  finding  the  smallest  value  of  s(p)  such  that 
y»p  has  order  p.  This  minimal  value  is  given  by 
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(5.2.3) 


s(p)  - 


p - 1,  2, 3,  4 


p + 1 p - 5,  6 
p + 2 p - 7 

unknown  p 2 8 


For  msthods  of  order  greater  than  seven,  a gap  develops.  For  instance,  eighth-order 
methods  with  eleven  stages  exist,  and  it  is  known  that  any  eighth-order  method 
requires  at  least  ten  stages.  For  arbitrary  p 2 8,  the  best  bounds  known  for  the 


optimum  value  of  s(p)  are 


(5.2.4) 


p ♦ d(p)  £ s(p)  S (p2  - 7p  + 14)  / 2 , 


where  d(p)  2 c In  p for  all  p sufficiently  large  (for  some  c > 0).  The  lower  bound  is 
given  in  Butcher  [75J  the  proof  is  quite  involved,  and  the  result  is  not  much  bette- 
than  the  "trivial"  lower  bound  s(p)  2 p (Hindmarsh  [74],  page  84).  A class  *cyRK  °* 
methods  such  that  *p  requires  only  (p2  - 7p  ♦ 14)  / 2 stages  is  given  in  Cooper  and 
Verner  [72]. 

We  first  consider  lower  bounds  on  the  complexity  C(p,or)  using  LRK  methods. 
The  "trivial"  lower  bound  s(p)  2 p will  be  used,  since  the  term  $(p)  will  be  small  when 
p is  small  and  will  not  affect  the  asymptotic  behavior  of  optimal  order  and  complexity 
for  p large.  It  i*  K-own  (Butcher  [64])  that  at  loast  0(p2)  of  the  subdiagonal  elements 
of  the  matrix  A (whose  elements  are  the  X]j  in  (5.2.2))  must  be  non-zero  in  order  for  A 
to  define  a p^-order  method.  Thus  there  exists  > 0 such  that 


(5.25) 


since  s(p)  2 p,  we  see  that 
(5.2.6) 


where  we  now  write 


d(p)  2 aL  pz  i 


e(82p(v»  2 N eL(v)  p 


eL(v)  min  e(vj)  . 


*-^J 
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Thus  (5.2£)  and  (5.2.6)  show  that  a lower  bound  on  the  cost  per  step  for  *p  is  given 
by 

(5.2.7)  cl(p)  - aL  p2  + N ej_(v)  p . 

Theorem  5.2,1; 

C|_(p,«)  - [aL  p2  + N eL(v)  p]  e*/p  . 

Proof;  This  follows  immediately  from  (4.6)  and  (5.2.7).  | 

It  is  clear  that  ^(p)  :■  Ml  ^l  p2  ♦ N bl(v)  p]  e*/p  satisfies  (2.10)  and  (2.11). 
We  claim  th*»t  ^ yields  a Gl  satisfying  (2.15).  Indeed,  write 

fL(p)  - fj(p)f2(p)> 

where 

ft(p)  Ml  aL  p 

and 

fo(p)  s*  p + » , where  r :■  N 6l(v)  / »l  • 

Clearly  fj  yields  a Gj  satisfying  (2.15).  Since  f2 's  * l>n**r  polynomial  with  a negative 
zero,  it  may  be  shown  thr.t  f 2 yields  a G2  satisfying  (2.15).  Tnus  fL  yields  a Gl 
satisfying  (2.15b  in  fact,  we  have 

(5.2.8)  Gl(p)  - cjtp)  + G2(p)  - p [1  ♦ (1  ♦ ap’1)"1]  • 


This  leads  us  to 

Theorem  5.2.2: 

Cl#(«)  **  [Ml  8l  s2  / 4]  a2  as  a t co  . 

Proof;  From  (5.2.8),  we  see  that  Gl(p)  ~ 2 p as  p t oo.  Since  (2.10),  (2.11),  and 
(2.15)  hold,  Pl*(«)  tends  to  infinity  with  a.  Thus 

a ■ GL(pL*(o))  ~ 2 Pl*(«)  as  a t oo  , 
i.e.,  Pl#(«)  ~ er/2  as  a t oo.  The  result  now  follows  from  Theorem  5.2.1.  g 
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We  now  turn  to  upper  bounds  on  complexity.  The  class  $cvrj<  derived  *n 
Cooper  and  Vernei  [72]  has  two  deficiencies,  the  first  of  which  is  that  no  uniform 
upper  bound  on  *L(j(pptH)  ‘s  known  f°r  ^CVRK»  'n  addition,  the  combinatory  cost  for 
this  class  of  methods  is  Ofp^)  as  p T oo.  Instead,  we  turn  to  the  basic  sequence  *qrk 
discussed  in  Appendix  A.  There,  we  prove  that  there  is  an  Mu>0  such  that 

(5.2.9)  »G<i»p,h)  S (My  In  (p  + e)  h)P  , 

provided  h £ hp,  where  hp  - 0((ln  p)“^)  as  p t oo.  Furthermore,  there  are  a large 
number  of  extra  zeros  in  the  matrix  A for  ?p  « #cRK*  *^s'n8  the  notation  of  Appendix 
A,  we  see  that  the  number  of  non-zero  entries  in  A is 

rs  t.  . tP~^  i2  + □ 

zi-0«i  2i-i  1 + P 

- p3/3  - p2/2  + 7p/6 

S p3/3  + 2p2/3 

for  p € Z ++.  Finally,  note  that  the  number  of  stages  s(p)  required  for  ?p  € #cRK 's 

(5.2.10)  s(p)  - [(p2  - 2p  + 4)/2J  S p2/2  + p 

for  p € Z ++,  which  shows  that  the  number  of  stages  required  for  a p^-order  method 
in  *CRK  asymptotically  equals  the  number  requires  for  a p^-order  method  in  *qvRK- 
Thus  (considering  the  combinatory  costs),  the  class  $CVRK  actually  costs  more  per 
step  than  does  #CRK;  'Sn0r’ng  the  combinatory  costs  would  have  caused  us  to  reach 
the  opposite  conclusion. 

First,  we  look  at  the  cost  per  step.  By  (5.2.10), we  see  that 

1 o 

(5.2.1 1)  e(5JUv)}  S - (p2  ♦ p)  N e,  Kv) , 

p 2 


where 


ey(v)  :«  max  i^isN  e(vj)  • 


Since  we  are  using  $crk>  *t  is  easy  t°  see  that  there  is  a bj  Z 2/3  such  that 

(5.2.12)  d(p)  5 (p3/3  ♦ bj  p2)  • 2N  . 
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Combining  (5.2.11)  and  (5.2.12),,  we  sea  that  the  total  combinatory  cost  per  step  Is 
bounded  by 

(5.2.13)  cy(p)  «*  N [ 2p3/3  + p2  ♦ $2  P 1 « 

where 

fil  :-  ey(v)  / 2 + 2 by  and  <S2  «l/v>  / 2 • 

Using  (5.2.9)  and  (5.2.13)  gives 
Theorem  5.2.3: 

Cy(p,a)  - IA.  N [2p3/3  * p2  + 02  p J In  (p  + e)  o*/P  . Q 
Now  we  look  at  the  optimality  theory  for  tne  upper  bound. 

Theorem  5.2.4.: 

(1.)  For  all  * > 0,  there  exists  py*{«)  such  that  (3.17)  holds. 

(2.)  py*(«)  increases  monotonically  with  e,  and 
Py*(«)  ~ a/3  as  e T co  . 

(3.)  Cy#(«)  increases  monotonically  with  a,  and 

Cy*(«)  ~ [ 2 My  N e3  / 31  ] a3  In  « as  a T co. 

(4.)  hy*(«)  - ( My  e3  In  o )"*  as  a T co  . 

Proof:  We  write 


fpjtp)  tm.  My  In  (p  + e)  cy(p) 

in  the  form 

fyfp)  - fi<P)f2(p)» 

where 

f j(p)  “ My  N p In  (p  + e)  and  f2(p)  » 2p2/3  + p + 02  • 

As  was  pointed  out  in  Section  2,  fj  satisfies  the  hypotheses  of  Theorem  2.1.  Now  we 
consider  f2>  Clearly  f2  has  no  positive  zeros?  it  may  be  seen  that  the  condition 
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bU*  2/3  implies  that  f2  has  a positive  discriminant  and  hence  has  no  complex  roots. 
Thus  t'l  has  only  negative  roots;  one  may  the>v  show  that  this  guarantees  that  f£ 
satisfies  the  hypotheses  of  Theorem  2.1.  By  Lemma  2.1,  the  same  may  be  said  for  f. 

Thus  py*  and  Cy*  behave  as  described  in  (1.)  of  Theorem  3.3.  We  also  see  that 
Gy(p)  ~ 3 p as  p t co.  Thus  the  estimate  in  (2.)  holds,  from  which  we  get  the  estimates 
in  (3.)  and  (4.).  | 

So  in  the  class  of  linear  Runge-Xuit?  methods,  we  find  that 

<5.2.14)  CL*(«>  - Of.2)  £ C*(«)  S Cy*(«>  - CX«3  In  c) 

as  a tends  to  infinity;  hence,  the  ratio 

Cy*(«)/C L*<«)  - Qim  In  «) 

indicates  the  gap  in  our  Knowledge  of  the  complexity  of  linear  Runge-Kutta  methods. 

Finally,  we  wish  to  compare  the  classes  of  Taylor  series  methods  and  LRK 
methods.  Write  Cy^y*  , C^y*  , and  Cy*  (respectively,  » ^l,LRK*  » anc^  Cy^*) 

for  Cy*,  Cy*,  and  C*  in  the  class  *y  (respectively,  the  class  Since  we  have 

only  asymptotic  expressions  for  these  quantities,  we  are  forced  to  use  an  asymptotic 
comparison.  If  f,  g : R++  -+  R**  satisfy  lim  aToo  f(a)  *•  lim  g(a)  - -Ho,  we  will 
write 


(5.2.15)  f < g iff  f(o)  « o(g(<r))  as  a t co  ; 

we  say  that  f is  asymptotically  less  than  g.  If  f < g,  there  is  an  > 0 such  that 
f(ct)  < g(a)  for  a > «q,  so  there  is  a non-asymptotic  interpretation  of  the  order 


T 


relation  <.  In  addition,  we  see  that  if  f < g,  then  g(«)  grows  much  more  quickly  than 
f(a)  does  as  a increases.  Using  the  results  of  (5.1.20)  and  (5.2.14),  we  then  have  the 
following 


Theorem  5.2.5:  Suppose  that  (5.1.16)  holds. 


(U 

If  N ■ 1,  then 

CU,T*  < CUIRK*  * 

(2.) 

If  N - 2,  then 

^T*  < culrk*  • 

(3.) 

If  N - 3,  then 

Cu,T*(«)  “ ^C^y^fe)) 
end 


Cu,LRK^*>  “ 

as  a T oo. 

(4.)  If  N £ 4,  then  Cy^y^*  ^ ^L,T*  * S 
If  (5.1.16)  does  02l  hold,  then  (1.),  (2.),  and  (3.)  may  be  false,  but  <4)  v-'H! 
certainly  be  true.  As  an  immediate  corollary  to  the  above  theorem,  we  have 
Theorem  5.2.6: 

(1.)  If  N « 1 and  (5.1.16)  holds,  then  CT*  < 0 y^*  . 

(2.)  if  N «s  4,  then  Cy^*  ^ Cj*  . | 

So  if  de-vatives  are  cheap  to  evaluate,  we  see  that  the  best  Taylor  series 
method  known  is  better  than  the  best  linear  Runge-Kutta  method  possible  for  the 
scalar  case  N - 1?  but  if  N £ 4,  the  best  linear  Runge-Kutta  method  known  is  better 
than  the  best  Taylor  series  method  possible. 


I 


Section  6 


Nonlinear  Runge-Kutta  Methods  for  the  Scalar  Case 

We  now  consider  a generalization  of  the  familiar  LRK  methods,  based  on  the 
description  In  Brent  [74],  [76}  A basic  sequence  said  to  be  a class  of 

nonlinear  ftunge-Kutta  methods  (abbreviated,  "NW<  methods")  if  each  increment 
function  may  be  written  in  the  form 

(6.1)  *p(X}»h)  rs(X[,h;  kg,  „ ^s-1) , 

where 

(6.2)  Kj  v(yj),  yj  Tj(x,-,h;  kg,  - ,kj_l)  (0  5 j S s - 1) 

for  suitable  functions  tj:  R^Rx(R^)t  -»  R^  (0  S j S sb  as  in  the  linear  case,  s - s(p) 
is  the  number  of  stages  (i.e,  evaluations  of  v)  of  *p.  Again,  for  noiational  convenience, 
we  do  not  indicate  the  dependence  of  the  Kj,  yj,  and  rj  on  Xj  and  p. 

In  the  remainder  of  this  section,  we  will  only  consider  the  scalar  case  N - 1, 
since  it  is  not  known  whether  NRK  methods  exist  for  larger  yalues  of  PI  In  this  case, 
(5.1.5)  shows  how  an  s-stage  NRX  method  of  order  p may  derived  from  a (p+li^-order 
iterative  method  for  solving  the  nonlinear  equation 

(6.3)  F(z)  - 0, 
using  Brent-Information  (Meersman  [76])  of  the  form 

(6.4)  $lB JF)  (F^F^F^^-.F^)}  . 

Brent  [741  [75]  used  this  transformation  to  derive  a sequence  of  (modified) 


Brent-Runze-Kutta  method?  ("BRK  methods"),  in  which  the  s-stege  method  has  orde' 
(65)  p - 2s  - 1 . 


an 


Furthermore,  Meersman  [“8]  proved  that  this  order  is  the  greatest  possible  in  the 
class  of  ail  such  BRK  methods.  We  now  extend  Meersman’s  result  to  include  all  NRK 
methods. 

Theorem  6.1:  No  s-stage  NRK  method  can  have  order  greater  than  2s  - 1. 

Proof:  Let  be  an  s-stage  method  with  order  p.  We  will  construct  (from  *)  an 
iterative  method  ^ of  order  q :■  p + 1 for  finding  a simple  zero  f of  an  arbitrary 
analytic  function  F : R -»  R. 

The  method  is  defined  as  follows.  Let  xg  be  an  approximation  to  f such  that 
F'  !s  nonzero  between  xg  and  f*.  (Since  F'(f)  j*  0,  such  an  xg  exists.)  Write  to  :■* 
F(xq);  without  loss  of  generality,  assume  tQ  < 0.  Now  apply  one  step  of  using  a 
step-size  -tQ,  to  the  problem 

xii)  - F'(xtt))-1  (t0  < t < 0)  with  x(t0)  - x0  , 
whose  solution  is  the  functional  inverse  nf  F 

x(o>  - rh o)  - r ; 

then  ^ is  given  by 

f(x0)  x0  - t0  ?(x0,-tg)  • 

By  the  defin',:on  of  order  for  iterative  methods,  it  is  clear  that  f has  order  q; 
moreover,  ^ uses  the  generalized  Brent  information  (Definition  IL3.3  of  Meersman  [76]) 

$GB,s  =-  {F(x0),F'(y0),F'(y1),..,F'(ys.1)}  . 

Suppose  that  yQ  r1  xq;  then  q i 2s  by  Theorem  IL3.3  of  Meersman  [76].  On  the  other 
hand,  if  yg  « xg,  then  ^ uses  the  Brent-information  (6.4);  by  Theorem  IL2.4  of 
Meersman  [76]  (also  due  to  WofniaKowsKi),  we  have  q £ 2s  in  this  case  also.  Thus  in 
either  case,  we  find  that 

p + 1 - q S 2s  , 


and  the  desired  result  follows.  B 
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.i ! 


Thus  ♦j^jgRK  is  informationally-optimal  in  the  class  of  NRK  methods,  in  the  sense 
that  each  fp  in  uses  the  minimum  number  of  stages  possible  for  a p^-order 

NRK  method. 

We  will  now  derive  lower  bounds  for  the  complexity  C(p,«)  via  NRK  methods. 
Clearly,  Theorem  6.1  implies  that 

(6.6)  e<S2p(v))  £ e(v)  (p  + 1)  / 2 , 
and  a linear  lower  bound  on  the  combinatory  cost  states  that 

(6.7)  d{p)  2 a^p 

for  some  a^  > 0.  By  (6.6)  and  (6.7),  a lower  bound  on  the  cost  per  step  for  *p  is 


(6.8) 


cL(p)  - (aL  + e(v)/2)  p + e(v)/2  , 


i 


which  leads  to 

Theorem  6,2:  CL(p,«)  - Mj_  [(aL  + e(v)/2)  p + e(v)/2]  e*/P  . 

Proof:  This  follows  immediately  from  (4.6)  and  (6.8).  | 

Note  that  fy(p)  :«•  Mj_  cy(p)  is  a linear  polynomial  with  a negative  zero;  it  then 
folio* i that  fj_  satisfies  the  conditions  of  Theorem  3.2.  So,  the  optimality  theory  of 
Sec.‘:o  , 3 holds;  in  particular,  we  have 

Theorem  6.3:  C|_*(a)  ~ e [a^  + e(v)/2]  a as  o t co  . 

Proof:  From  (4.7)j  and  (6.8),  we  find  that 

GL(p)  - p2  / (p  + 0_i) , where  0 :■  1+2  aj  / e(v) , 
and  so  Gl(p)  ~ p as  p t ooj  thus  Pl*(«)  * « as  a T co.  The  result  follows  by  letting 
p - P(_*(a)  in  the  definition  of  C^p^r).  | 

Next,  we  consider  upper  bounds  on  the  number  of  operations  required.  Instead 
of  using  *MBRK>  w®  use  c*ass  *BRK  °f  "unmodified"  BRK  methods  described  in 
Appendix  B,  where  it  is  shown  that  there  is  an  My  > 0 such  that 


(6.3) 


<rG(*p.h)  5 (MU  h>P  i 
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no  such  bound  is  Known  for  ♦|w$0r«.  In  addition,  *mbrk  requires  the  solution  of  p - 1 
linear  systems  of  equations,  the  i^  having  p - i unknowns,  in  order  to  perform  a 
“reorthogonaiization."  So  the  smallest  Known  combinatory  cost  for  this  class  is  about 
0(p3-8I)  arithmetic  operations}  this  is  obtained  by  using  Strassen’s  technique  for  linear 
systems  (described  in  Borodin  and  Munro  [75]}.  On  the  other  hand,  most  of  the 
combinatory  cost  for  in  is  involved  in  finding  the  coefficients  of  the 

polynomial  pn+j  (see  Appendix  Bh  once  these  coefficients  are  Known,  the  remaining 
combinatory  cost  is  0(p  In  p)  as  p T oo.  An  estimate  of  how  much  work  is  required  to 
compute  these  coefficients  is  given  in 

Iggma  giis  Let  xg,  yj,  _ , yr,  w0,  z0,  _ , zr  be  given,  and  let 

Qfx)  l!^1  q,  x* 

be  the  unique  polynomial  of  degree  at  most  r ♦ 1 satisfying 

CKx0>  - w0  , Q'(x0)  - z0 , and  Q'ty,)  - Zj  (1  S i s r)  . 

If  T(r)  is  the  time  required  to  compute  qg,  , qr+j,  then 

T(r)  - 0(r  In^r)  as  r t oo  . 

Proof:  The  coefficients  qj,  2q£,  , (r+Dqr+^  of  Q'  may  be  computed  in  time 

0(r  In^r)  by  using  a fast  algorithm  for  computing  the  coefficients  of  the  Lagrange 
polynomial  interpolating  the  points  (xq,  z q),  (y^^),  ~ , (yr,zr)j  see  Borodin  and  Munro 
[75]  for  details.  Then  Q(r)  ope.  ations  yield  qj,  ...  , qr+j^,  and  Horner’s  rule  gives  qg 
with  0(r)  additional  operations.  | 

Thus  there  exists  ay  > 0 such  that 

(6.10)  d(p)  S ay  p ln^(p+e)  . 

In  order  to  simplify  matters  a bit,  note  that  Theorem  B.1  implies  that 


(6.11) 


e(Slp(v))  S e(v)  p . 
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Although  the  estimate  above  is  not  exnct  for  p > 2,  it  is  asymptotically  equal  to  that  in 
Theorem  B.l.  (If  necessary,  the  sharper  estimate  given  there  may  be  used,  but  the 
calculation  of  optimal  order  involves  considerably  more  detail,  the  results  of  which  are 
not  particularly  enlightening.)  Combining  (6.10)  and  (6.11),  we  see  that  the  cost  per 
step  is  bounded  by 

(6.12)  cu<p)  - e(v)  p + ay  p ln2(p+e)  . 

Thus  (6.9)  and  (6.12)  imply 

Theorem  6.4:  Cy{p,cr)  - My  [e(v>  p + ay  p ln2(p+e)]  e“/p  . | 

We  now  determine  the  behavior  of  Py*(«).  Here  we  find  that  fy(p)  :■  My  cy(p) 
may  be  decomposed  as 

fy(p)  - fi(p)f2(p), 

where 

f l(p)  My  e(v)  p and  f2(p)  1 + 0 ln2(p+e) , 

and  0 ay  / e(v)  . Clearly  fj  and  f2  satisfy  (2.10)  and  (2.11),  and  f]^  yields  a Gj 
satisfying  (2.15).  We  need  only  check  that  f2  yields  a G2  satisfying  (2.15).  But 
G2(p)  - 2 0 p2  In  (p+e)  / D2(p) , where  D2(p)  (p+e)  f2(p) , 

so  that  setting 

g2(p)  :«  0 p ln2(p+e)  [in  (p+e)  - 1]  + 2 0 e ln2(p+e)  + (p  + 2e)  In  (p+e)  + p , 
we  find  that  p > 0 implies 

G2'(p)  - 2 0 p g2(p)  / D2(p)2  > 0 . 

Thus  the  optimality  theory  applies. 


Theorem  6.5: 

(1.)  py*(e)  increases  monotonicaily  with  e,  and 
Py*{«)  ~ a as  * t oo  . 

(2.)  Cy*{«)  increases  monotonicaily  with  «,  and 

Cy*(w)  **  My  ay  e * ln2a  as  * t oo  . 

(3.)  hy*(*>  - (My  e)-i  as  a t oo  . 

Proof;  (1.)  follows  from  the  fact  that  Gy(p)  ~ p as  p t oo;  (2.)  and  (3.)  follow 
from  (1.)  and  Theorem  6.4.  i 

So  in  the  class  of  nonlinear  Runge-Kutta  methods,  we  find  that 
(6.13)  CL*(«)  - 0(a)  6 C*(e)  £ Cy*{«)  - Of«  in2*) 

as  a tends  to  infinity;  so,  the  ratio 

Cy*(«)  / CL‘<*>  - CKIn2e)  as  * t oo 

indicates  the  gap  in  our  Knowledge  of  the  complexity  of  nonlinear  Runge-Kutta 
methods. 

Finally,  we  wish  to  compare  the  classes  of  Taylor  series  methods  and  NRK 
methods.  Adopting  tha  notation  at  the  end  of  Section  5.2  in  an  obvious  manner,  we 


Theorem  6.6:  If  (5.1.16)  holds,  then  Cyj*  ^U,NRK*  • 

Proof:  Immediate  from  (5.1.20)  and  (6.13).  | 

Thus  if  derivatives  of  v are  easy  to  evaluate,  the  best  Taylor  series  method 
Known  is  better  than  the  best  nonlinear  Runge-Kutta  method  Known.  However,  if  the 
cost  of  evaluating  the  K^  derivative  of  v increases  faster  than  0(ln  K)  as  K T oo,  then  it 
is  easy  to  show  that  the  opposite  will  be  true. 


Section  7 


Numerical  Results 


a 


In  the  previous  sections,  we  computed  (for  several  classes  of  methods)  that 
order  which  minimized  the  work  required  to  attain  a given  error  criterion.  Here,  we 
consider  actual  numerical  results  of  optimal  order  and  minimal  cost  for  various  test 
problems  and  classes  of  methods.  The  optimal  order  for  a given  error  criterion  was 
determined  by  finding,  for  each  method  implemented,  the  coarsest  mesh  that  allowed 
the  error  criterion  to  be  satisfied;  the  resulting  complexities  were  then  compared  to 
determine  the  optimal  order.  The  error  measure  used  was  the  "endpoint  error,"  i.e., 
the  co-norm  (see  e.g.,  Stewart  [73],  pg.  164)  of  the  difference  between  the  true  and 
computed  solutions,  evaluated  at  the  endpoint  of  the  interval  of  interest  (the  unit 
interval  I).  All  testing  was  carried  out  on  the  Carnegie-Mellon  University  Computer 
Science  Department’s  PDP-10  in  ALGOL  and  FORTRAN,  using  double  precision. 

The  first  problems  considered  were  of  the  form 

(7.1)  x(t)  - X x(t)  x(0)  - 1 

on  the  unit  interval  L Although  this  problem  is  easy  to  handle  analytically,  any  general 
problem  of  the  form  (2.1)  may  be  locally  approximated  by  a linear  system  of  ore  nary 
differential  equations  (see  e.g.,  Hindmarsh  [74],  pp.  1 7 - i 8).  If  tL»  coefficient  matrix  of 
this  linear  system  is  diagonalizable,  an  uncoupled  set  of  scalar  equations  of  the  form 
(7.1)  will  result. 

These  problems  were  solved  via  Taylor  series  methods;  the  optimal  order  is 
given  in  Table  7.1  for  the  choices  of  X indicated.  Here  the  optimal  order  was  taken  to 
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be  that  order  which  minimizes  the  number  of  evaluations  of  the  right-hand  side  of  (7.1) 
required  to  attain  the  desired  ermr  criterion.  As  expected,  the  number  of  evaluations 
required  increases  as  the  errt  'iterion  i decreases.  Moreover,  the  optimal  order 
also  increases  monotonically  as  « decreases,  just  as  the  theory  predicts. 

We  next  turn  to  the  solution  of  the  test  problem 
(7.2)  x(t)  - cos2x(t)  x(0)  - 0 . 

For  this  problem,  we  searched  for  the  optimum  "unmodified"  Brent-Runge-Kutta 
method.  For  this  problem,  the  optimal  order  was  taken  to  be  that  for  which  the  actual 

! 

j CPU  time  (in  milliseconds)  required  to  solve  the  problem  to  within  a given  s was 

minimized.  Since  there  is  a certain  amount  of  randomness  in  such  a measure,  the  mean 
time  for  ten  runs  was  analyzed.  Not  surprisingly,  it  turned  out  that  the  order  which 
minimized  the  CPU  time  also  minimized  the  number  of  evaluations  of  the  right-hand 
side  of  (7.2).  Since  the  (n  + 2),h-order  method  requires  the  zeros  of  the  Jacobi 
polynomial  Gn(2,  2,  • ),  and  the  best  set  of  values  available  only  contained  the  zeros 
for  1 £ n s 8 (Table  25.8  of  Abramowltz  and  Stegun  [64]),  only  the  methods  of  order 
not  exceeding  ten  were  implemented. 

The  results  for  problem  (7.2)  are  given  in  Table  7.2.  Here,  the  optimal  order  p*, 
the  optimal  number  of  mesh  points  n*,  the  minimal  number  of  evaluations  C*,  and  the 
minimal  mean  CPU  time  C*  are  given.  Note  that  these  all  behave  as  predicted.  In 
addition,  we  computed  tho  ratio  of  the  mean  CPU  time  for  a fourth-order  method 
C*(4,  • ) to  the  minimal  mean  runtime.  As  the  theory  predicts,  this  ratio  appears  to  be 

S 

t 

increasing  without  bound  as  i tends  to  zero.  (The  same  behavior  was  found  for  the 
ratio  Ce(4,  • ) / C*  , where  Ce(4,  • ) is  the  number  of  evaluations  required  by  a fourth- 
order  method.) 


Finally,  we  looked  at  the  "hard"  problem 


5s^»  wSft^sa  *= 
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(7.3) 


X|<t)  - XJ.J  orjj(t)  xj(t)  Xj(t)  (1  S I S 2) 

*jj(t)  - Y|j  J*“  exp(-i»jj  (t  - r)2)  r-1  sin  r dr  (1  £ i,  j £ 2) 
(where  "exp"  denotes  the  exponential  function),  with  initial  conditions 

xj^(O)  ■ X2(0)  - 1 

The  ?jj  were  all  taken  to  be  one,  while  tne  rjj  were  taken  to  be 


*11  » 1 , >12  “ * 22  * » **21  ” 10'^ 

(This  system  of  differential  equations  is  similar  to  the  system  governing  a two-species 
gas  chemical  reaction;  see  e.g.,  Finlayson  [71].) 

Since  the  system  (7.3)  is  nonscalar  and  nonautonomous,  the  Brent-Runge-Kutta 
methods  are  not  appropriate.  Since  the  derivatives  of  Xj(t)  are  not  readily  available, 
the  Taylor  series  methods  are  not  particularly  easy  to  apply.  Thus  we  used  linear 
Runge-Kutta  methods  for  the  solution  of  (7.3).  The  particular  methods  RKp  of  order  p 
(1  £ p £ 8)  used  were  as  follows. 


RK1  ...  Euler’s  method 

RK2  ...  Ralston  [66]  (5.6-40)  "modified  Euler" 

RK3  ...  Ralston  [66]  (5.6-45) 

RK4  . . . Ralston  [66]  (5.6-48)  "classical  method" 

RK5  . . . Cassity  [66] 

RK6  . . . Butcher  [64b]  (first  method  on  page  192} 

RK7  . . Shanks  [66] 

RK8  . . . Cooper  and  Verner  [72] 

The  methods  of  order  less  than  eight  have  the  optimal  number  of  stages  per  step, 
while  the  method  of  Cooper  and  Verner  has  the  minimum  number  of  stages  of  all 
eighth-order  methods  known  (see  Section  5.2). 

Most  of  the  work  involved  in  solving  (7.3)  was  in  evaluating  ajj(t).  An  obvious 
change  of  variable  reduces  this  to  a Gauss-Hermite  quadrature;  a twonty-point 
quadrature  (Table  25.10  of  Abramowitz  and  Stegun  [64])  was  used  for  maximal 
accuracy.  The  Chebyshev  rational  function  approximation  given  on  page  356  of 
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FrS'oerg  [69]  was  used  to  compute  (sin  r)  / r for  |r|  £ 1;  the  system  doubie-precision 
sine  routine  was  used  for  jf  j > 1 . 

Since  so  much  of  the  time  required  to  solve  (7.3)  was  spent  in  evaluating  «jj(t), 
the  measure  of  cost  was  the  number  of  evaluations  of  the  set  {ffjj(t) : 1 £ I,  j £ 2}  j 
that  is,  we  measured  the  number  of  evaluations  of  the  (vector)  right-hand  side  of  (7.3). 
(Moreover,  the  amount  of  computer  time  required  to  search  for  the  optimum  was  so 
great  ns  to  preclude  running  the  problem  a large  number  of  times  and  averaging  the 
results,  as  was  done  in  the  previous  example.)  Results  are  given  in  Table  7.3,  where 
p*,  n*,  and  C*  (defined  as  for  (7.2))  are  given  as  a function  of  the  error  criterion.  The 
table  stops  at  t - 10"^,  since  the  eighth-order  method  (i.e.,  the  highest-order  method 
implemented  for  testing)  was  reached  at  that  level.  Again,  note  that  the  theoretical 
results  predicted  are  confirmed  in  this  difficult  examp.e. 

So,  our  three  numerical  examples  yield  data  which  agree  with  the  theoretical 
result  that  the  optimal  order  p*(a)  increases  with  a - In  « . Moreover,  in  Sections  5 
and  6,  we  saw  that  p*(<r)  - 0{«)  as  e t oo?  i.en  the  optimal  order  increases  linearly 
with  a.  The  data  in  Tables  7. 1-7.3  support  this  result. 

Further  testing  still  remains  to  be  done.  In  these  examples,  we  picked  problems 
that  were  well-suited  to  one  particular  type  of  method  (e.g^  it  was  easy  to  get  the 
derivative  information  required  by  Taylor  series  methods  for  (7.1)).  Future  testing 
should  look  at  problems  that  are  "neutral"  in  the  sense  that  the  informations  required 
for  the  various  classes  of  methods  are  equally  hard  to  obtain.  This  would  allow  the 
comparison  of  various  classes  of  methods.  In  addition,  we  point  out  that  "fast" 
methods  of  polynomial  manipulation  wer„  not  used  (due  to  the  additional  programming 
involved  in  designing  such  a package);  perhaps  such  a package  should  be  Implemented 
for  future  testing  of  the  Taylor  series  and  nonlinear  Runge-Kutta  methods. 
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TABLE  7.1 

Taylor  Series  Methods  for  Test  Problem 
x(t)  - X x(t)  x(0)  « 1 


-Iogi0i 

X ••  -e 

X ■ -1 

X - -1/e 

X ■ 1/e 

X » 1 

X - e 

1 

2 

3 

1 

1 

3 

8 

2 

9 

4 

2 

2 

4 

9 

3 

11 

G 

3 

3 

G 

11 

4 

12 

7 

4 

4 

7 

12 

5 

14 

8 

5 

5 

8 

14 

6 

15 

9 

G 

G 

9 

15 

7 

1G 

10 

7 

7 

ie 

1G 

8 

17 

11 

8 

8 

11 

18 

9 

19 

12 

9 

9 

12 

19 

Notes: 

1.  In  all  cases  except  X - -e,  i « 10'  *,  the  optimal  mesh-size  was 
h • 1.0;  for  this  exceptional  case,  it  was  h - 0.5  . 

2.  Entry  in  table  is  the  optimal  order  for  the  given  X and  i.  This  equals 
the  minimal  number  of  function  evaluations  required  to  solve  the 
problem  on  the  entire  unit  interval,  except  for  the  exceptional  case 
noted  above,  where  four  was  the  minimal  number  of  evaluations. 


TABLE  7.3 

Linear  Runge-Kutta  Methods  for  Test  Problem 
x(t)  « J*al  aij(t)  xjtl)  xj(t)  Xj(0)  - 1 (1  5 » 5 2) 

) “ Yjj  exp{-rj:  (t  - r)2)  t'1  sin  t dr  (15  i,  j 5 2) 


1 
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C* 
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8 
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Section  8 

Summary  and  Conclusions 


In  this  thesis,  W9  have  constructed  s methodology  for  studying  the 
computational  complexity  of  one-step  method*  for  the  numerical  solution  of  ordinary 
differentia!  equation  initial-value  problems.  We  devolved  lower  and  upper  bound*  on 
the  complexity  of  « given  msthod,  and  showed  how  to  pick  that  method  within  a basic 
sequence  which  rsirslmtee*  complexity.  Under  very  general  hypotheses  {which  were 
later  verified  for  a number  of  commonly-used  classes  of  methods),  we  saw  that  the 
optimal  order  increases  as  the  error  criterion  t decreases,  tending  to  infinity  as  t 
tends  to  zero;  this  is  in  contrast  to  the  situation  in  iterative  complexity  (Traub  and 
WoSnsskowski  [76]).  Moreover,  in  many  cf  the  specific  classes  of  methods  studied,  we 
saw  that  the  optimal  step-size  does  just,  tend  to  zero  witn  s,  a result  indicating  that  it 
is  important  to  not  assume  that  h tends  to  zero.  These  result  of  optimal  order  and 
stsp-$iz&  were  then  used  to  find  bounds  on  the  complexity  of  tiolving  the  equation  to 
within  z given  «,  using  a gtv%n  class  of  methods;  using  these  bounds,  we  were  able  to 
compare  the  “goodness"  of  several  such  classes. 

Vfe  now  turn  to  some  issues  that  have  been  r-'rso  by  this  study.  Probably  the 
most  important  point  is  that  we  have  found  evidsn.s  that  high-orier  methods  may  be 
of  practical  (i.e,  computational),  ss  well  ss  theoretical,  hierest.  However,  we  need  to 
learn  much  more  about  them,  the  crucial  point  not  being  that  of  gatting  maximal  order 
for  a given  information  set,  but  that  of  getting  minimal  complexity  for  a method  of 
given  order.  For  example,  the  optimally-ordered  dess  ♦mbrx  has  greater  complexity 
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than  the  class  indeed,  the  Differences  in  combinatory  cost  far  outweigh  any 

advantages  the  former  has  over  the  latter,  as  was  pointed  out  in  Section  6.  (Even  if 
we  were  to  take  the  drastic  step  of  ignoring  combinatory  cost,  we  would  still  be  faced 
with  the  fact  that  is  known  to  be  order-convergent,  while  no  such  result  is 

known  for  4ygR|<-)  A similar  situation  arises  when  we  compare  the  classes  and 

4>CVRK  of  linear  Rungo-Kutta  methods  (Section  5.2). 

Finally,  we  consider  some  open  questions. 

(1.)  In  a number  of  instances,  we  have  only  been  able  to  show  "trivial"  lower 
bounds,  i.e.,  lower  bounds  which  are  linear  in  the  amount  of  information  needed. 
However,  the  best  algorithms  known  have  complexity  which  grows  faster  than  linearly 
in  the  size  of  the  information  set.  How  may  we  narrow  this  gap?  (Note  that  this  is  an 
issue  that  touches  almost  all  areas  of  complexity  theory.) 

(2.)  What  are  the  possibilities  of  extending  this  ana',  sis  to  include 
"polyalgorithms"  for  the  solution  of  initial-value  problems?  It  is  often  wise  to  vary  the 
step-size  from  stto  to  step;  furthermore,  many  existing  programs  allow  the  order  to 
vary.  It  would  be  useful  to  have  a complexity  theory  that  includes  these  methods. 

(3.)  We  assumed  throughout  this  thesis  that  infinite-precision  real  arithmetic 
was  available.  It  would  be  of  great  interest  to  study  the  complexity  of  initial-value 
problems  using  variable-precision  arithmetic,  a far  more  realistic  model. 

(4.)  What  is  the  complexity  of  using  multistep  methods  to  solve  initial-value 
problems?  We  have  some  preliminary  results  for  the  class  of  Adams-Bashforth 
predictor/Adams-Moulton  corrector  methods,  using  Sgpx  *°  f'n<*  the  necessary  starting 
values;  this  work  will  be  reported  in  a future  paper.  (This  class  is  order-convergent, 
and  a result  similar  to  Theorem  3.3  holds;  it  also  appears  that  the  choice  of  starting 
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method  is  of  great  importance  in  the  complexity  analysis.)  Besides  the  multistep 
methods,  many  oth,  classes  of  methods  remain  to  be  analyzed,  such  as  extrapolation 
methods,  spline  methods,  multistep  Runge-Kutta  methods,  and  special  methods  for 
"stiff"  equations.  (Of  course  this  list  is  by  no  means  exhaustive;  see  Gear  [71]  or 
Hindmarsh  [.  4]  for  further  discussion.) 

(5.)  The  error  equation  (3.1),  (3,2)  holds  for  non-iterative  methods  for  the 
solution  of  a number  of  other  problems  such  as  boundary-value  problems  for  ordinary 
and  partial  differential  equations.  Hence,  we  suspect  that  a great  deal  of  the  analysis 
in  Sections  2,  3,  and  4 may  go  through  unchanged.  A long-term  goal  is  the  study  of 
such  problems  from  a complexity  viewpoint. 


Appendix  A 

Error  Bounds  for  a Basic  Sequence  of 
Cooper-Runge-Kutta  Methods 


In  this  Appendix,  we  describe  a subclass  of  a class  of  linear  Runge-Kutta  ("LRK") 
methods  due  to  Cooper  [69}  We  shall  first  prove  the  following 

Theorem  A.1:  There  is  a basic  sequence  $crk/  ^-RK  methods  such  that 
(1.)  Each  ?p  € *crk'  re9u'res 

s(p)  (p2  - p + 2)  / 2 
evaluations  of  v per  step. 

(2.)  There  exists  an  My  > 0 such  that 

(A.1)  'G<Vh)  5 (Hi  In  (p+«)  h)P 

for  h < hp  - 0((ln  p)"*). 

We  use  the  notation  of  Cooper  and  Verner  [72}  Let  p < Z be  given;  define 


p:  Z + n [0,  p]  -♦  Z + by 


(A.2)  p(j) 


2k«0  K B ^ 2 if  i * P 
s if  j - P . 


where  we  write  "s"  for  "s(p)"  as  defined  above.  Next,  a set  {£q  , _.  , of  integers  is 
defined  by  picking  $g  p,  and  setting  £j  (i  + 0)  to  be  the  unique  integer  in  [1,  p] 


satisfying 


p(fe  - 1)  < i S p(fc) 


We  now  pick  Uq  , ...  , us  ( I satisfying 


Uq  - 0 , u$  - 1 , Uj  i1  0 if  i i1  0 


*V1 


s? 


t 


} 

i 


K 


! 


. '-.  fc£Jf  V7fc  A 
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(A.5)  ({j  - {j  «nd  i i<  j)  implies  Uj  f*  Uj  . 

Finally,  we  pick  a matrix  of  coefficients  A :«  {fy:  0 £ j 5 i-1,  1 S I £ s}  such  that 
(A.6)  X,j  = 0 if  ' $j  < |j  - 1 (1  S i,  j S s> 

and 

(A.7)  X,-|  uf  - (r+ir^Ujr+1  (0  s t s $ - 1, 1 <£  i S s)  . 

Cooper  ano  Verner  [72]  point  out  that  these  conditions  may  always  be  fulfilled;  the 
resulting  A defines  a p**1 -order  LRK  method  with  s stages. 

We  are  interested  in  a choice  of  uq  , _ , u,  which  will  give  a small  error 
coefficient.  To  this  end,  we  will  choose 

(A.8)  {uj:  (j  - n}  - {(1  + xkn)  / 2 : 1 5 k S n)  (lSr<Sp-l), 
where  x^n  , ...  , xnn  are  the  zeros  of  the  Jacobi  polynomial  Pn  Pn^»^  (see 
Szegb  [59]>.  Since  these  zeros  are  distinct  and  lie  in  [-1, 1] , conditions  (A.4)  and  (A.5) 
may  be  satisfied. 

Now  we  are  able  to  exhibit  a solution  to  the  i^  system  in  (A.7).  First,  note  that 
4.e  equation  for  r - 0 may  be  separated  from  the  others,  since  Uq  - 0.  Setting 

n :■  |j  " 1 , 

we  see  that 

(A.9)  Xj0  - Uj  - - uj  - I { X|j:  j < i and  Z n } , 

the  last  by  (A.6).  We  wish  to  determine  the  nonzero  Xjj,  i.e.,  those  Xjj  for  which  £j  2:  n 
and  j < i.  So  setting 

Xjj  -0  unless  j<  {j[  , ~ , jn) , 
we  see  that  the  remaining  Xjj  are  the  solution  of  the  system 

2k-l  V \ " (r+irl  uir*1  (1  * r 5 n)  * 

Thus  the  Xjj  are  the  weights  for  an  interpolatory  quadrature  formula  on  [0,  Uj]  with 


thal  0 < <>kn  s */2.  Usin8  (8.9.2)  of  SzegB  [59],  we  then  find 

Aikn  “ 0(^^n  JJ.  n+^  [Pn(cos  d)  / (cos  »J  - cos  tfkn)]  sin  $ d«>  . 

22S®.  L-  ^l>n+l  * ^i,n+l  s ^k,n+l/2*  We  consider  the  integral  over 
[^ln/2»  ^i.n+l] » since  Theorem  15.4  of  SzegB  [59]  proves  that 
0<k5/2n-3)  [|  J'  | . | J*l"/2  |)  . 0(„-l,  . 

here  tr.e  integrand  is  the  same  as  in  the  preceding  integrJ.  But  the  proof  of  (15.4.12) 
in  SzegB  [59]  extends  almost  immediately  to  a proof  that  the  remaining  integral  is 
0(k~  n),  since  (15.4.12)  is  proved  by  order-of-magnitude  estimates.  Thus  pjkn  - 
0(n“l)  ■ 0(n"*  In  n)  for  Case  1. 

~ —■  - ^k,n+l/2  * ^i,n+l  s 3^k,n+l^2*  We  consider  the  integral  over 

t^kn/2’  ^i.n+l^  * since  SzegB  [59]  shows  that 

0(k5/2n-3)  | f*  | - 0(n_1)  . 

^kn'2 

As  in  (15.4.13)  of  SzegB  [591  w®  have 

jjj$  • 0(nk-3/2)  I,  ♦ l2  . 

Here 

with 

D(d)  [cos  (Nt>  -T  y)  - cos  (Nl>kn  ♦ <y)]  / [cos  I?  - COS  dkn]  , 
where  N :■  n ♦ 3/2  and  y :•»  -3f/4,  and 


'2  - S*$  Rn<».'W  S'n  * " ' 0(nk'3/^, ' 

with  Rn  the  remainder  term  in  (8.8.2)  of  SzegB  [59].  Unfortunately,  the  proof  that 
(15.4.14)  of  SzegB  [59]  is  bounded  does  not,  extend  to  a proof  that  is  bounded, 
since  the  proof  of  the  former  requires  that  the  interval  of  integration  be  symmetric 
about  0kn.  However,  it  is  straightforward  to  verify  that 

- (XI)  |sin  Ni>  / $|  dd  - 0 (In  n)  . 

Thus  Mjkn  ■ CXn“2k  In  n)  « 0(n-*  In  n)  for  Case  2. 

Cc.se  2-  3tfk,n+l  5 ^i.n+l  5 3»/4.  «*e  consider  the  integral  over 

[3<Jkn/2,  0 j n+j] , since  SzegB  [59]  proves  that 

°<k5/2n'3'  I fonn/2  I - • 

But  the  proof  of  (15.4.19)  in  SzegB  [59]  extends  to  prove  that  the  refining  integral  is 
0(k‘5/2n)  (as  in  Case  1).  Thus  nikn  • CXn-1)  - CKn^  In  n)  for  Case  3. 

Case  4:  3r/4  5 $in+1  £ *„+l,n+l*  W#  considar  the  inte8ral  over 

[3w/4,  i>in+1]  , since  SzegB  [59]  shows  that 

0(k5/2n-3)  | Jgx/4  j - Oin”1)  . 

As  in  Cases  l and  3,  the  proof  of  the  above  may  be  extended  to  pfi.ve  a similar  bound 
on  the  integral  of  interest.  Thus  >ijkn  ■ 0(n”^)  * 0(n  ^ In  n)  in  Case  4,  completing  the 
proof  of  the  Lemma  g 

Thus  (A.9)  and  Lemma  A1  show  the  existence  of  a X > 0 such  that 


(A.10) 


IXjjl  S X In  (fj  ♦ e) ; 


here  X is  independent  of  p.  Moreover,  tha  result  for  the  case  i » s may  be  sharpened. 
We  see  that  X^  i 0,  since  the  Uj  for  the  s*h  system  in  (A.7)  are  the  abscissae  for 
Lobatto  quadrature.  Thus 


(A.l  1) 


iSji  ■ $ ».i  • >• 


the  consistency  condition  in  the  last  equality  being  a consequence  of  (A.7)  with  r ■ 0. 
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Proof  of  Theorem  A.1:  As  in  Cooper  and  Verner  [72],  we  define 

*j  x(ujh)  - Kj 
and 

«i  •-  Jo1  «"W  du  - >ij  i<Ujh) 

for  0 £ i £ s;  note  that  ig  ■ ig  “ 0.  Let  z(h)  be  the  computed  approximation  to  x(h)j 


(A.  12) 


h*1  ||x(h)  - z(h)||  - ||  h’1  [x(h>  - x(0)]  - 2?^  Xj,  kj  || 

s iiyi  * in-Jo  h,  -ill 


S ||6S||  + max  „ p-1  ||ij|| , 

the  last  by  (A.6)  and  (A.11).  By  the  analyticity  of  x,  there  is  an  Aj  > 0 such  that 
Uj  s-  h'1  ||  x(ujh)  - Z*Lg  (Ujh)rx^)(0)  / r!  ||  £ (Ax  h)*1 
and 

tjj  II  x(ujh)  - Z^o1  (Ujh f x<*>0)  / r!  ||  S h)*1  , 
so  that  the  definition  of  5j  gives 


(A.  13) 


ll«ill  ^ 0j  + 2™  |Xjj|  tij 

S (At  h)5i  + iXjjl  (Ax  h)*1 
S (Agh)*' 


for  a suitable  A2  > 0.  Thus  (A.1 2)  becomes 

(A.1  A)  h_i  Ux(h)  - z(h)||  s (A2  h)P  + max  £.  „ p-1  ||i,||  . 

We  now  use  Lemma  l.l  of  Cooper  and  Verner  [72]  and  (A.6)  to  find  that  if  L is  a 
Lipschitz  constant  for  v,  then  there  exists  A3  > 0 such  that 
||ij||  S hL  116,11  ♦ hL  Z':l0  IXjjl  max  j ||(j|| 

S (A3  h)^'  1 + (A3  h)  In  (£j  e)  max  j ||i||| , 

the  last  by  (A.10)  and  (A.13)i  here,  the  maximum  is  taken  over  all  j < i such  that 
£ £j  - 1.  A straightforward  induction  shows  that  if  (1  + In  2)  A3  h < 1,  then 


\ 
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ll*ill  * <A4  In  (Jj  + e)  h)fi+1 
for  a suitable  A4  > 0.  Combining  this  with  (A.14),  w©  find 

<A.15)  h"1  ||x(h)  - z(h)|{  S (A5  In  (p+«)  h>P  , 

the  desired  bound  for  the  local  error  for  a single  unit  step. 

To  extend  (A.15)  to  a global  error  result,  we  must  look  at  the  Upschitz  constants 
for  the  increment  functions.  Let  L be  a bound  on  ||Vv||,  and  write  "Wp(y»h)"  to 
indicate  gradient  with  respect  to  the  vector  variable  y.  Now 

||V„p<y,h)||  S Ip}  |y  max  ||7kj(y,h)|| 

- max  o$i£s-l  HVkj(y,h)|| , 

where  we  write  "kj(y,h)"  to  indicate  the  dependence  of  kj  upon  y and  h.  By  the 
definition  of  kj(y,h),  we  find 

Vkj(y,h)  - W(u>  [lNxN  + h Xjj  Vkj(y,h)] , 
where  u y + h ij^  Xjj  kj(y,h)  and  lfvjKN  '8  *n  N*N  '*ntity  matrix.  Taking  norms  in 
the  above  gives  the  result 

f i S LX  + hLX  [ In  (£,+©)  max  { fj:  j < i and  f j fc  (j  - 1 } ] , 
where  fj  :«  ||Vkj(y,h)||.  Writing  Xp  for  the  Upschitz  constant  for  *p,  it  is  easy  to  see 
that  (A.16)  and  the  above  inequality  imply 

Xp  5 Ip}  (hLX)*  np*  In  (p+e-k) , 

which  is  bounded  for  all  p,  provided  that  h S hp  < (LX  In  (p+e))"*.  Thuc  (A.1)  follow* 
from  this  result,  (A.15),  and  Theorem  3.3  of  Henrici  [62].  | 

The  value  for  s(p)  indicated  in  Theorem  A.1  may  be  improved  somewhat  by 
noting  that  since  we  are  using  a Lobatto  quadrature,  higher  order  may  be  expected 
with  fewer  steps.  Indeed,  if  we  use  the  strategy  outlined  in  the  comments  following 
Theorem  4 of  Cooper  and  Verner  [72],  we  have 
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Theorem  Aj2:  There  exists  a basic  sequence  #CRK  °*  l-RK  methods  such  that 
(A.l)  holds  and  requires 

s(p)  :-  L(p2  - 2p  + 4)  / 2J 


evaluations  of  v per  step.  8 


Ip  WWIW  W*  BHWfTO  *' ,f  l»  WOTt 


t, 


Appendix  B 

Order-Convergence  of  a Basic  Sequence  of 
Brent-Runge-Kutta  Methods 


In  this  Appendix,  wo  describe  a subclass  of  a class  of  iterative  methods  for  the 
solution  of  scalar  nonlinear  equations.  This  subclass  will  then  be  used  to  generate  an 
order -convergent  basic  sequence  of  nonlinear  Runge-Kutta  methods. 

Lemma  EJL:  Let  F:  DcR  -»  R have  a simple  zero  f,  and  suppose  that  F is 
analytic  at  f.  Pick  k,  m f ZH  with  m + 1 i k.  Then  there  is  a sequence 
*km  !“  *^kmn  : n € Z4"*"}  of  stationary  multipoint  methods  without  memory  such  that 
th*?  following  hold: 

(1.)  The  method  uses  the  information 

«kmn  {F<x0),~,F<'n)(x0),F<k)(yi),-.,F<k)(yn)} 

(the  points  yj  , _ , yn  being  suitably  chosen)  to  compute  a new 
approximation  xj  to  f from  a given  approximation  xq  by  setting 

X1  *kmn<xQ>  • 

(2.)  There  exists  aB>Oar>danhg>0  such  that  if  |xq  - f|  S h^  , then 
|xj  - f|  S (B  |x0  - f|>*  for  all  n < Z^, 

where 

(8.1)  p min  (m  ♦ 2n  ♦ i,  2m  ♦ n ♦ 1)  . 

Before  proving  the  Lemma,  we  describe  how  the  method  a'^mn  computes  an 
improved  approximation  xj  from  the  old  approximation  xq  . 


Algorithm  for  computing  Xt  ^Kmn^‘ 

(1.)  Let8  |F(x0)/F'(x0)| . 

(2.)  Let  Zj  be  an  approximate  zero  of 

Pl<x)  X^q  (x  - x0)'  F«>(x0)  / i! 

satisfying 

(B.2)  zj  » x0  ♦ 0<i)  and  !p1(z1)(  S (Aj  8)m+1  , 
where  Aj  is  independent  of  n. 


(3.)  Let 


yj  x0  + cjn(z«  x0)  (1  S i S n) , 


where 


r M 


“in  ^ + xif«'  I ^ 

and  x^n  > ...  > xnn  are  the  zeros  of  the  Jacobi  polynomial 
Pn(x)  p(h-l,  m+l-K)(x) 

(see  Szegb  [59]), 

(4.)  Let  pn+j  be  the  polynomial  of  degree  at  most  m + n that  interpolates  the 
information  n , and  let  xj  be  an  approximate  zero  of  pn+^  satisfying 
(B.3)  xj  - xq  ♦ 0(5 ) and  lpn+l^xl^  s (A2  • 

where  A2  is  independent  of  n and  p is  given  by  (B.1). 

Here  we  use  the  notation  of  Brent  [74].  Clearly,  < C'(k,  m,  n),  the  only 
difference  being  that  conditions  (B.2)  and  (B.3)  replace  (2.2)  and  (2.4)  of  Brent  [74].  It 
is  easy  tc  soe  that  (B.?j  and  (B.3)  may  be  reaped  by  using  flog2(m+l)]  - 1 and 
riog2(p/(m+l))]  iterations  of  Newton’s  method,  with  the  respective  starting 
approximations  of  xq  - F(xq)  / F'(xq)  and  Zj  . 

Proof  qf  Lemma  B.l:  Lot  xj'  be  the  exact  zero  of  pn+j  near  xq.  We  then  find 
that  there  is  a { between  x^'  and  z^  such  that 
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(8.4)  |F(xjO|  S IPn+l<2i>  “ f<21>1  * iPn+1^  " t“l'  ' zll  * 

Using  (6.3),  the  analytidty  of  F,  and  standard  techniques  of  Interpolation  theory 
(Tr  8ub  [64]},  it  is  easy  to  show  that  (2.S)  end  (2.10)  of  Brent  [74]  may  fee  rewritten  es 

ipn+i(x)  - Rx)|  S (Agl)^1  and 

(B.5)  m+_ 

|p'+1(x)  - F'(x>|  S (A4Dm+" 

for  |x  - Xq|  S 41.  (Hare  all  constants  A,,  wlil  be  independent  of  n.)  Slm'.lariy,  we  find 
that 

|xi'  - f|  £ (A5  l)m+fl  and  iz^  - f|  S (A$  <>m+l  * 
so  that  the  triangle  inequality  gives 

(B.6)  jx.'  -z{\  S (A7l)m+l  • 


Using  (B.4),  (B.5),  and  (£6),  we  see  that 


|f(V)l  * lpn+l<2i>  ' F<21>1  * <AS  5)2m+n+l 

6 |pn+i(Zi)  - Fi<21>I  * lF2^2l^  + *A8  *)2m+n+l  • 


v'here 


Fj(x)  (x  - ^(xq)  / l!  arKl  F2{x>  F(x)  " Fl(x> 

Clearly  1F2(x>1  S (A9  |)m+2n+l  . so  that  (0.7)  becomes 

(as)  |F(xj#/i  s lpn+ifci>*Fi(zi,l  + <A10I)'  • 

As  in  Brent  [74],  we  now  write 

pn+l(x)  • ri(x)  + r2^  1 

where  r,  (i  - l,  2)  is  the  polynomial  of  degree  at  most  m n satisfying 

r,W(x0)  - F,0>(x0)  (0  < j i m) 


If  we  let 


r,<k)(y;)  - Fj(K)(yj)  (i  S j S n)  . 


P(x)  rj(x  ♦ xq)  - FA(x  ♦ x9) » 
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and  write  t :■  zj  - xq  {in  this  Appendix  only),  we  find  that 

P(i)(0)  * 0 (0  <;  i S m)  and  P(k)(o-m  .)  - 0 (1  £ i S n)  . 

We  may  easily  alter  the  proof  of  Lemma  4.3  in  Brent  [74]  to  show  that 

rl<zl>  " Fl<zl>  - P<«>  - 0 . 

Thus  (B.8)  becomes 

(B.9)  IRxj')!  < |r2<zi)l  ♦ (A10  5y»- 

To  bound  the  remaining  term,  let  us  write 

'Zl>»  ■ . 

recalling  that  r2  has  a zero  of  multiplicity  m at  xq  . Using  the  notation  of  Stewart  [731 
we  see  that  the  nonzero  coefficients  of  r 2 are  given  by  the  solution  of  the  linear 
system 

Wy  - c , 

where 

«jj  a!'1  (1  £ i,  j £ n) , 

:«  aj+fn  «i+m  (j  ♦ m):  / (j  + m - K)!  (1  £ j £ n) , and 

T,  *k  F2(k)(y,)  / a^"k+1  (MiSn)  . 

SifVie  WT  is  a Vandermonde  matrix,  we  find  that  the  entries  of  U * W~*  are  given  by 

**  ' rn-i,n-l,j  ^ ^r?<j  ^®jn  ' ®rn^ » 

where 

^n-l.j  2 “pj.n  - aPM,n  • 

the  sum  being  taken  over  all  multi-indices  pj  ...  p^  not  including  j (Gregory  and 
Karney  [69]).  Since  there  are  fewer  than  2n  summands,  each  of  which  lies  in  [0,  1]  , 
we  see  that  £ 2n  , implying  that 

|vjj|  £ 2 / nrMj  («jn  - «rn)  . 


So  we  have 


<aio> 


where 


111  ■ ®h  Vi' 

S n 2"  m»  lsisn  | .k  F2W(y,)  / [«j»*  <V<«jn>l  I , 


G^x)  G^m  ♦ 1,  in  ♦ 2 - k,  x)  - E"^  (x  - «,.„) 


(see  Abramowitz  and  Stsgun  [64]}. 


Now  it  is  dear  that 


, . m-K  , . m-k 
m*M  Sj<nW«jri  - l/«nn 


By  Theorem  8.9.1  of  Szegb  [591  may  show  that 

«nn  * An  n-2; 

using  this  result  and  (22.5,2)  of  Abramowitz  and  Stegun  [641  we  find  that 


mfiX  lSjln  VVr 


0.11) 


S A12n2<m^  ^ + 1J  max  lsj£n  |Pn'<xjn)i'1  . 


By  the  symmetry  relation  (4,1.3)  of  Szegfi  [591  w®  may  assume  that  0 S xjn  < 1.  Using 
Theorem  8.9.1  of  Szeg&  [591  we  may  show  that 

lfV<*jn>r1  S (Ai3>n. 

and  so  (B.10),  (B.  11),  the  definition  of  F2,  and  the  above  imply  that 

fcjl  S (Al4  «)m+2r‘*1  , 


yielding  the  result 


|r2(Zi)|  S S"  aj+m  *i+m  S n max  lstSn  fc,i  S (Alg  8)' 


m+2n+l 


So  (B.9)  becomes 


F(*iO|  S (Al6  6)* 


By  Taylor's  Theorem,  this  implies 


ixj'-ri  * <a17  iy»  . 

The  desired  result  then  follows  from  (6.3)  and  from  (25)  of  Srent  [741  8 


'"t'’ 


*3*^^fi*  V. 


We  now  describe  the  basic  sequence  • The  methods  in  this  basic  sequence 


are  given  by 


pjtxQ^h)  v(xq), 

1*2^0  > h)  v(xq  + h t *)  i 


and  for  p £ 2, 


*p<*o  • h)  h_1  W'i,i,p-2(xo>  - xo3  * 

with  l,p-2  aPP*'ed  th®  function  F given  by  (5.1.5}  and  the  approximation  x^  to 
xj^  being  given  oy  an  appropriate  number  of  iterations  of  Newton’s  method  (as 
described  above). 

Theorem  BJL  The  basic  sequence  *brk  is  order-convergent  with  respect  to  the 
global  error.  Moreover,  the  number  of  stages  s(p)  required  by  ( $qrk  is  8‘ven  by 


s(p)  - 


p if  p S 2 

p - i if  p > 2 . 


Proof:  We  use  the  notation  of  Lemma  8.1,  writing  z(h)  for  the  computed  p^- 
order  approximation  x^  to  x(h)  and  pn+j(  * , xq)  for  the  polynomial  pn+i  . The  result 


of  Lemma  B.1  is  that 


tT1  |z(h)  - x(h)I  S (Bh)P  , 


the  desired  result  for  a single  unit  step.  To  prove  the  global  result,  we  must  consider 
the  Lipschitz  constants  for  $gRK. 

We  implicitly  differentiate  the  result  pn+j(x^',  xq)  ■ 0 to  find 


where 


?p<xo  • h)  " -b’1  Qn+l(xl'-  x0)  * ‘p^O* 


Qn*l<xl'.  x0}  " 1 * d2  Pr.*l<xr*  x0>  / dl  Pn+l(xl'>  x0> 


a 


•p^o*  " h"  <d/dxo>  ■ xi/l 


t 

t 

} 

1 

t 

I 


i 
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It  is  easy  to  see  that  and  x^'  are  analytic  functions  of  xq  . Sine®  their  difference 
tends  to  zero  uniformly  on  the  domain  of  v as  p t co  , it  follows  that 

ptoo  ‘p(*0^  " ® • 

We  claim  that 

Qn+i^i',  Xq)  - Oth  Inn)  as  n T oo  , 

uniformly  in  xq  . Tc  see  this,  note  that  we  may  write  the  interpolation  polynomial  Pn+i 
in  terms  ''f  Jacobi  polynomial  Pn  , finding  that 

Pn4l<^0>  - (-Dn  (h/2)  Pn(t)  dt  4 h 'Kx0)lJ-1Ikn  - h, 

where 

ffx)  2 (x  - xq)  / [h  v(xq)]  - 1 

and 

'kn  •-  [2  <1  ‘ *kn>  v<Vk)  Pn^knlf1  « * » *V»  / « - *kn>  d»  • 

Nr  , 

*1  Pn4]<xl'.x0>  - HJHPnttpMxo)  ♦ (i  * fj)  2 g(*Kn>  , 

where 

fi  s-  rtxjo , 

Lkn(x)  FnM  / [Pn'(>W  (x  * xkn>3  * ™d 

g(t)  1 / [U  4 f)  v(xQ  4 (1  4 t)  h v(x0)  / 2)]  . 

By  (8.21.10)  of  SzegS  [591  the  first  term  in  the  expression  for  dj  Pn+i(xj'  , xq)  goes 
to  zero  as  n T oo  . A minor  modification  of  the  proof  of  Theorem  14.4  of  Szegfi  [59] 
show-;  that  the  s jr.i  in  the  remaining  term  tends  to  g(f(x(h)))  as  n t oo  . So 

dl  pn4l(xi' » x0^  ~ vWh))-1  as  n t oo  . 

Using  Lemma  A.l  and  techniques  similar  to  ihose  yielding  the  above  estimate,  we,  find 
<^2  Pp+itxi' , Xq)  - 0(h  In  n)  - v(x(h))'^  as  n T oo  . 

This  gives  the  estimate  ciaimed  for  Qn+jtxj' , xq)  . 


So  the  Lipschitz  constant  for  < ♦qrk  grows  as  the  logarithm  of  p.  By 
Proposition  4.3,  #brk is  order-convergent.  | 
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